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Co3qaH B Weax HHPOPMUpOBaHHA YNTATeJIbCKOM ayAUTOPHU O HOBeEMINMX JOCTWKCHHAX HM MepcieKTHBaXx B OOaCTH 
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Based on Three-Dimensional Equations of Elasticity Theory 
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Abstract 

Introduction. Functionally graded materials are of great use, because heterogeneity of properties enables to control the 
strength and rigidity of structures. This has caused great interest in the topic in the world scientific literature. The 
construction of solutions to such problems depends significantly on the type of boundary conditions. In this paper, we 
consider the equilibrium of a thin-walled circular cylinder whose mechanical properties change along the radius. 
Homogeneous boundary conditions were set on cylindrical surfaces that had not been considered before, the effect was 
on the ends. The mathematical formulation of the problem was carried out in the linear theory of elasticity in the 
framework of axisymmetric deformation. Expressions were constructed for the components of the stress-strain state of 
the cylinder, in which some coefficients were found from the solution to the resulting system of linear algebraic equations. 
Materials and Methods. The material of the cylinder was linearly elastic, the elastic modulus of which depended linearly 
on the radial coordinate. The basic research method was the asymptotic method, in which half the logarithm of the ratio 
of the outer and inner radii acted as a small parameter. Iterative processes were used to construct the characteristics of the 
stress-strain state of the cylinder. 

Results. Homogeneous solutions to the boundary value problem were obtained for a linearly elastic functionally gradient 
hollow thin-walled cylinder. An analysis of these solutions made it possible to reveal the nature of the stress-strain state 
in the cylinder wall. For this purpose, an asymptotic analysis of the solutions was carried out, relations for displacements 
and stresses were obtained. It was determined that those solutions corresponded to the boundary layer, while their first 
terms determined Saint-Venant edge effect similar to the plate theory. 

Discussion and Conclusion. The analytical solution to the equilibrium problem of a thin-walled cylinder inhomogeneous 
in radius constructed by asymptotic expansion can be used for numerical solution to a specific problem. For this, it is 
required to solve the obtained systems of linear algebraic equations and determine the corresponding coefficients. The 
resulting asymptotic representations provide analyzing the three-dimensional stress-strain state. The selection of the 
number of expansion terms makes it possible to calculate displacements and stresses with a given degree of accuracy. 


This analysis can be useful in assessing the adequacy of applied calculation methods used in engineering practice. 


Keywords: linear theory of elasticity, functionally graded material, thin-walled hollow cylinder, homogeneous solutions, 


boundary layer, variational principle 
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AHHOTalna 

Beedenue, DyHKMoHaIbHO-rpamqMeHTHble MaTepHasIbl HAXOLAT OOMbIIOe IPHMeHeHHe, T.K. HEOMHOPOAHOCTh CBOMCTB 
HO3BOACT YIPaBIATh MPOUHOCThIO MW %KCCTKOCTHIO KOHTPYKIMH. ITHM BbI3BAH OOIbIIOM HHTepec K WaHHOM Teme B 
MUpoOBOM Hay4dHoH WuTepatype. locrpoeHue pellieHua TakuX 3ajad CyNIeCTBCHHO 3aBHCHT OT THIa TpaHH4HBIX yCJIOBHI. 
B Hactosileti padoTe paccmaTpHBaeTca paBHOBeCHe TOHKOCTeHHOrO KpyrOBOrO WWIMHApa, MexaHWueckHe CBOMCTBa 
KOTOPOTO 34MCHAIOTCA BAONb paguyca. Ha WuMHApH4yecKkHXx MOBePXHOCTAX 3aaHbI OMHOPOAHbIe TPaHWGHble YCOBUA, 
KOTOpbIe WO ITOTO He paccMaTpHBAaJINCh, BO3elCTBHe OKa3bIBaeTCA Ha Topulax. MarematTw4eckaad NocTaHoBKa 3aa4qu 
OcyllecTBIAeTCA B JMHeMHOM TeopHu yipyrocTH B paMKax OCeCHMMeTpHYHON Aedopmaiuuu. B padote mocTpoeHsI 
BbIP@KEHHA JI KOMIOHCHT HalipsKeHHO-eOpMMpoBaHHOTO COCTOAHUA IWIMHApa, B KOTOPBIX HeKOTOpbIe 
KOSPPUIUMCHTHI HAXOAATCA M3 PeMICHHA MOWYICHHOM CHCTeMBI JIMHeMHBIX alreOpav4eckux ypaBHeHHii. 

Mamepuaaoi u memooot. Martepual WunuHypa ABIIAeTCA JMHeEMHO yIPyruM, MOLyJIb ypyrocTu KoTOporo JIMHeiHHO 
3aBHCHT OT payvaIbHOM KOOpAMHAaTbI. OCHOBHBIM MeTOJOM HCCIeOBaHHA ABJIAeCTCA ACCHMITOTH4eCKHH MeTOL, B 
KOTOPOM B KayecTBe MaJIOrO MapaMeTpa BbICTyHaeT MOMOBHHA JOrapHdMa OTHOINeCHHA BHeINHerO UH BHYTpeHHero 
panuycos. Jia NocTpoeHua XapakTepHcTHK Hallps2%KeHHO-eOpMMpoBaHHOrO COCTOAHUA IMIMHApa TMpHMeHeHbI 
UTepallMOHHBle IIpOleccsl. 

Pezynomamet ucciedosanua. Jina sMHewHO-yupyroro (dyHKIMOHaIbHO-rpagqWeHTHOTO MONOrO TOHKOCTeHHOTO 
IWIN Apa WOUYYeCHEI OTHOPOAHbIe pellieHHa KpaeBoli 3ayaun. AHasIH3 3TUX PeIICHuM MO3BOIMeT PacKPbITh XapakTep 
Hallps.KeHHO-edOpMHpOBaHHOTO COCTOAHHA B CTeHKe IMIMHApa. C ITO Lebo MpoBeseH ACHMITOTHYeCKHH aHasH3 
pelleHuii, MOJyUeHbI COOTHOMICHHA JIA MepeMelleHHi MU HallpsyKeHHM. YcTaHOBJICHO, UTO ITH PeIIICHHA COOTBETCTBYIOT 
MOrpaHHYHOMy CJIOIO, Ip ITOM HX MepBble YWIeHbI OMpeeuAOT KpaeBont 9dext Cen-Benana, aHasOrM4Hblli Teopuu 
TMT. 

O6écyotcoenue u 3axu04enue. TloctpoeHHoe Cc MOMOLIbIO ACHMIITOTHYECKOLO pa3zIOXKeHHA aAHAIMTMYeCKOe pelleHue 
3alauH O paBHOBeCHH HeOHOpOAHOrO 0 paswycy TOHKOCTeHHOrO IMIMHEpa MO#KeT ObITb HCIONb30BaHO IIA 
YMCJICHHOrO pelieHHA KOHKpeTHOM 3aqauu. Jina 9TOrO HYyKHO pelIHTb TMOyYeHHbIe CHCTeMBI JIMHeMHBIX 
anreOpawuecKux ypaBHeHHii WM OlpeseuMTb coorBeTcTByroulMe Kosddunnents. TlomyaeHHple acuMmToTH4eckHe 
IIpeJicTaBICHHA TO3BOUAIOT aHaIM3HPOBaTb TPeXMepHOe HalipsxKeHHO-ed@opMupoBaHHoe cocTosHHe. Bsroop 
KOJIM4YECTBA YICHOB pa3zIO2XKCHUA MO3BOIAeCT PACCUHTATb WepeMelCHHA H HalIpMKCHHA C 3a/JaHHOM CTeMeHbIO TOUHOCTH. 
OTOT aHaiv3 MOxKeT OBITh MONe3eH MPH OWeHKe aleKBAaTHOCTH UpPWKaqHbIX MeTOOB pacueTa, IPHMeCHAeMBIX B 


MHKeHeEpHOH lpakTuKe. 


KsoyeBbie CioBa: JMHelHadt Teopua yupyrocTu, PyHKUHOHANIbHO-rpaqMeHTHbIM MaTepHasl, TOHKOCTeHHbIM MOJIbIM 


UWHJIMVHAp, OMHOPOAHbIe pelwleHHusA, nMorpaHvuHEit con, BapHaliMOHHbIit TIPHHUAIT 


BaarogxapHocrnu: aBTOPbI BbIpaxKaroT OaroapHOcTh pedakWHu U peleH3enHTy 3a BHHMAaTeJIBHOe OTHOLMICHHE K CTaTbe UH 


TIPCAIO7KCHHA, KOTOPbIe MO3BOJIMJIM HOBbICHTb Ce KavyecTBO. 


Ismayilova J. Analysis of Stress-Strain State of a Cylinder with Variable Elasticity 
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Introduction. Functionally graded materials are widely used in various constructions. Due to the dependence of 
mechanical properties on coordinates, it is possible to control the stress-strain state (SSS) of parts. An example of using 
such inhomogeneity is a cylinder whose mechanical properties depend on the radius. In this case, the cylinder may be of 
interest as a separate structure, or as being a division subring of a compound body, e.g., connecting two media with widely 
different properties. When calculating the SSS of a thin-walled cylinder, some applied theories can be used, an assessment 
of their adequacy. Especially in the case of inhomogeneous properties, it can be carried out using computer modeling or 
asymptotic analysis based on a three-dimensional formulation. The latter determines the relevance of this study. 

A number of studies have been devoted to the investigation of the SSS of hollow cylindrical bodies within the 
framework of the linear theory of elasticity. In [1, 2], the mechanical behavior of a radially inhomogeneous cylinder in a 
three-dimensional formulation was studied on the basis of the spline collocation method and the finite element method. 
In [3], the SSS of a cylinder whose properties depended on the radius loaded with uniform internal pressure was described. 
In [4], an analytical study was carried out for a functionally graded piezoelectric cylinder. In [5], an exact solution was 
constructed to a radially inhomogeneous hollow cylinder with exponential Young's modulus, with constant Poisson ratio 
and power Young's modulus. In [6, 7], an analytical solution to the axisymmetric thermoelasticity problem for a 
continuous cylinder was obtained using the direct integration method when the coefficient of linear thermal expansion 
was an arbitrary function of the radius. In [8], a general asymptotic theory of a transversally isotropic homogeneous 
hollow cylinder was developed. New groups of solutions were obtained for a transversally isotropic homogeneous 
cylinder. The comparison of the constructed solutions to the solutions constructed using applied calculation methods was 
given. In [9, 10], some boundary value problems of elasticity theory were studied for a functionally gradient isotropic and 
transversally isotropic (the plane of isotropy was perpendicular to the axis) cylinder, in the case when the elastic modules 
were arbitrary continuous functions of the radius of the cylinder. In [11], an analysis of the bending deformation problem 
for a radially inhomogeneous cylinder was carried out. The analysis of the above papers shows that not all types of 
boundary conditions on cylindrical surfaces have asymptotic representations of solutions. 

In this article, on the basis of an asymptotic analysis of three-dimensional equations of elasticity theory, the features 
of the SSS of a thin-walled cylinder whose properties vary linearly along the radius were studied. In this case, the inner 
border was fixed in the axial direction and was free in the radial direction. 

This goal was achieved using several steps: asymptotic integration of differential equations and the construction of 
homogeneous solutions; derivation of formulas for the components of the displacement vector and stress tensor; 
consideration of boundary conditions on the face surfaces. 


Materials and Methods. Radially inhomogeneous hollow thin-walled cylinder 
T= {r E iG ; Fal be [0; 2r], ze [=i Ly } is considered in a cylindrical coordinate system with the origin on its axis. The 
problem of its equilibrium, in the case of fixing cylindrical surfaces along the axis and zero normal stresses, is solved in 


an axisymmetric formulation under the action of stresses at its faces. 


The boundary value problem consists of the equilibrium equations [8]: 


Oo 0o 19 oe O45 


rT p 'Z 4 — 0) ‘ 1 

or Oz r ) 
09, i 09. re oO, 0, (2) 
or Oz r 


where 6,,,0,,. O4)>O,, — components of the stress tensor. 


rr? 


Defining relations [8]: 


s, = (2642) vi +) (3) 
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0, =(26+2) aH vi rte) (4) 
7 Oz or r 
u Ou. Ou 
6,, =(2G+A)— +A “4+ |, a) 
” ( ) r (= Oz a 
S. -6{& aa (6) 
. oz Or 


Here, uv, =u,(7,z), u, =u,(r,Z) — components of the displacement vector. 
Lamé parameters vary linearly along the radius: 
G(r) =G.,r, ArH dr, (7) 
where G,,A, — constants. 


After substituting (3)-(7) into equations (1), (2), the dimensionless system of equations takes the form: 


O2u Ou Ou, Ou, O2u 
(2G, +,)| —? +e— |+e(G, +A, Jew +1,62e? —+G,62e% —* —2G,6u, =0, (8) 
op? Op Opes 0§ 08? 
O2u Ou Oru Ou (ol 
G,| —= +e— |+(2G, +A, )e?| e% —++e" — |+(G, +A, )ee? —*2 =0. (9) 
op? Op 0S? 0§ Opog 
Here: 
1 r Zz ; ‘ : 1 Lh ; ‘ 
p=—lIn a new dimensionless coordinates; ¢ = a — | —in the case of wall thinness, small parameter; 
é Nt % q 
ii, Uu,. u. el G.r, : 
m=Vnn, pe[-bl], 6e[-Gl], /=2;u,=+, w=, Y= , G, =—*; G, —some parameter having the 
% % % G, G, 


dimension of stress. 


Let us consider a problem in which homogeneous boundary conditions are set on the lateral surfaces of the cylinder: 


Melo =O (10) 
Spl, =0- (1) 
Stresses are applied to the faces of the cylinder: 
Orefe uy = 4s (P)> (12) 
Oes|. 4, =Ls(P)» (13) 
(St). 
oO or O,. = Pr: O,. = ox dimensionless stresses 
pp G, 7 PE G, 2s . 


Stress vector components 1,,(p),t,,(P), (s=1;2) satisfy the equilibrium conditions. 


To construct homogeneous solutions, we seek the components of the displacement vector in the form: 


u,(psE)=u(p)ews, u,(p:5)=w(p)er. (14) 

Substituting representations (14) into the system (8)—(11), we obtain: 
(2G, +A, )(u"(p) +eu' (p)) + eae ((G, +A, )w'(p) +A ew(p)) + €2G, (a2e2 —2)u(p) =0, (15) 
G,(w"(p)+ew' (p))+(2G, +A, )e? (aeu(p)+ a2e2w(p))+e(G, +A, )aeru' (p) =0, (16) 
wo, = 9, (17) 
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=0. 


p=+1 


[ (2G, +o )u'(p)+ Ero (u (p)+ aew(p)) 


(18) 


We investigate boundary value problems (15)-(18) at se +0. To solve (15)-(18) at ¢ +0, we use the asymptotic 


method [9-13]. 


Nonzero solutions (15)—(18) correspond to the third iterative process, the components of the displacement vector are 


searched for in the form of expansions over a small parameter: 


i() Selec (peace) a), 
w®) (p) = (i (p) + ew, (p)+---), 
a=e!(B, +eB, +---). 


After substituting expansions (19) into equations (15)—(18) for terms of the first order, we have: 
(2G, a5 Ay U5 (P) +By (G, +o) Way (P) +G, aU39 (p) = 0, 
GW (p) +B, (G, + Do) Miho (p) + (2G, +o )B3 Wo (p) =0, 


=0, 


p=tl 


((2G, +A )usy (p) +A BoW5o (p)) 


Following [13], the spectral problem (20)—(23) corresponds to a potential solution for the plate. 


Wo (p) 


p=+1 . 


Thus, the solutions are presented in the form: 


a) UG) (p38) = exT, (-2G,B3, sinBy, sn(up)*+0(@))exp{ *(B +eBi, +e) s| . 


un (p38) - EDT, (2G, dx Si Bo, £05(Bp) +O(6))<exp{ * (Bu +B, + es) 6] : 


Here, B,, is the solution to the equation: 
cos Bo, =0. 


The stresses corresponding to solutions (24), (25) have the form: 


a” - DT, (—4G2B3, sin Box £05(Bup) +0(6))<exp{ * (Bu +B, +++) s] ’ 


os” = 3 T, (4G38;, sin By, sin(yp)+0(6))ex0| “(Bo +e€B,, + oo) | > 
kal € 


= 


3:1 < 300 1 
Gy y= 2 [, (4G28;, sin Bo, £05 (BP) +O(@)}exo{ “(Pa +eB,, +++) s| 
oi” - (6). 


b) uO) (p; é) - ey F, (2G,B3, cos By, cos (Bo:P) ss O(e)) “ exp( “(Pu +é€B,, ++ JB) 


i=l 


ue (p; 5) = ey F, (2G,B3, COS By, sin(Bo;P) + O(e)) “ exp( 4 (Bo, + eB, ++ ) s| , 


Here, B,, is the solution to the equation: 
sin B,; = 0. 


The stresses corresponding to solutions (31), (32) have the form: 


O82) = LF (—4G38}, cos By, sin (By,P) +O(c))es9{ “(By +éB,, +36} 


of = LF, (46283, cos B,, c05(Bup)+0(€))exP( 2 +eB,, +6} 


(19) 


(20) 
(21) 


(22) 


(23) 


(24) 


(25) 


(26) 


(27) 


(28) 


(29) 


(30) 


(31) 


(32) 


(33) 


(34) 


(35) 
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0° = 5 F (4638, cosB,, sin (Bp) +0(6))err( (Bs —— 8} ” 


of” =O(e). (37) 
The general solution (15)—(18) will be a superposition of solutions (24), (25), (31), (32): 


u, (p:&) =) 7, (-2G,Bi, sin Bo, sin (By,p) + O(e))x 
xexp( M(B +eB, + ++) g}robF (2G,B3, cos By; cos(By,P) + 0(@))exp( (Pn +éB,,+ oo) é} (38) 


Ue (p;6) = ey T, (2G,B%, sin, cos(By,P) + O(e))x 


> 
un 


«exp( 1B +B, +e) eS (2G,B3, cos By, sn(yp)+0(@) exp 2(By, +éB,, +98), (39) 


i=l 
Solutions (24), (25), (31), (32) have the character of a boundary layer. When moving away from the faces, solutions 


(24), (25), (31), (32) decrease exponentially. 
To determine constants 7,, F,, we use Lagrange variational principle. The variational principle takes the form [8]: 


i 


oy AL 
2 (6. fis OM, * (ox hs Jou, | Loe =0. (40) 
Substituting (24-36) into (40), we have: 
EM Tug = Paz: (41) 
$0,Fo = pl a 


Here: 


M x = 16G5B5; ok (i — Box i‘ sin Bo, sin Bo; sin (By, “Pao Bafa g{ Oa =Pul} (npu j# k) 
, € 
M ,, =16G;B3, sin? By, [exe ei 24) (npu j =k) 
1 l 
Po; = 2G}; sin Bo, it. (p)cos(B,,p)-t, (0)sin(B,p))do-evo{- a 
1 
Q;, = 16G5B5 Bo, (B,, — Boi y" COS By ; COS Bo, sin (By, —Bo: \e*[- Gest), Saal (npu i# j) 


2p, 2B, / 
Q, = 166585; cos’ By; feo{ 2). e( 74) (npu i=j ) 
, €& €& 


| oj! 
Poy) = 2G,B3; COSBy, ue (p)cos(B,,p) +12, (p)sin(B,,p))dp- cxn(- re }: 


+5 (tm (p)cos(B, ,P) —t,,(p)sin (B,,p))dp cx © 


Jj 
€ 


| . By jf 
+] fie (p)cos(B,,P) +t, (p)sin(B,,p))dp exe{"“}} 
TTF ely #4; 

F, = Fy +éF, Hees 


Ismayilova J. Analysis of Stress-Strain State of a Cylinder with Variable Elasticity 


Constants T,,,F,, (p =1, 2,...) are found from systems of linear algebraic equations (41), (42), analogous to which are 
studied in [13]. 

Research Results. Tp the article, in an axisymmetric formulation, the solution to the problem of linear elasticity theory 
for a functionally graded hollow thin-walled cylinder, whose properties vary in thickness according to a linear law, was 
considered. Homogeneous cross boundary conditions were set on the lateral surfaces of the cylinder, and a stress vector 
was set at the faces. The constructed homogeneous solutions satisfied boundary conditions on cylindrical surfaces. For 
their construction, an asymptotic approach based on the expansion by a small parameter characterizing the relative 
thickness of the cylinder was used. To account for inhomogeneous boundary conditions at the faces, systems of linear 
algebraic equations similar to those studied in the literature were obtained. It was shown that the constructed SSS solutions 
had a boundary-layer character, which corresponded to the edge effect similar to the theory of inhomogeneous plates, 
which bears the name of Saint-Venant. 

Discussion and Conclusion. Usually, when studying the SSS of thin-walled structures, applied calculation methods 
are built that reduce the dimension of the problem. In this regard, the task of determining the range of geometric and 
mechanical parameters in which these methods give acceptable accuracy is critical. The solutions to three-dimensional 
equations constructed in the work on the basis of asymptotic analysis make it possible to assess the adequacy of such 
applied theories with a predetermined accuracy threshold. In addition, these solutions can find application in the 


evaluation of numerical solutions to problems for structures with functionally graded materials. 
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Abstract 

Introduction. The widespread use of piezoelectric materials in various industries stimulates the study of their physical 
characteristics and determines the urgency of such research. In this case, modal analysis makes it possible to determine 
the operating frequency and the coefficient of electromechanical coupling of piezoelectric elements of various devices. 
These indicators are of serious theoretical and applied interest. The study was aimed at the development of numerical 
methods for solving the problem of determining resonance frequencies in a system of elastic bodies. To achieve this goal, 
we needed new approaches to the discretization of the problem based on the finite element method and the execution of 
the software implementation of the selected method in C# on the .net platform. Current solutions were created in the 
context of the ACELAN-COMPOS class library. The known methods of solving the generalized eigenvalue problem 
based on matrix inversion are not applicable to large-dimensional matrices. To overcome this limitation, the presented 
scientific work implemented the logic of constructing mass matrices and created software interfaces for exchanging data 
on eigenvalue problems with pre- and postprocessing modules. 

Materials and Methods. A platform was used to implement numerical methods .net and the C# programming language. 
Validation of the research results was carried out through comparing the values found with solutions obtained in well- 
known SAE packages (computer-aided engineering). The created routines were evaluated in terms of performance and 
applicability for large-scale tasks. Numerical experiments were carried out to validate new algorithms in small- 
dimensional problems that were solved by known methods in MATLAB. Next, the approach was tested on tasks with a 
large number of unknowns and taking into account the parallelization of individual operations. To avoid finding the 
inverse matrix, a modified Lanczos method was programmatically implemented. We examined the formats for storing 
matrices in RAM: triplets, CSR, CSC, Skyline. To solve a system of linear algebraic equations (SLAE), an iterative 
symmetric LQ method adapted to these storage formats was used. 

Results. New calculation modules integrated into the class library of the ACELAN-COMPOS complex were developed. 
Calculations were carried out to determine the applicability of various formats for storing sparse matrices in RAM and 
various methods for implementing operations with sparse matrices. The structure of stiffness matrices constructed for the 
same task, but with different renumbering of nodes of a finite element grid, was graphically visualized. In relation to the 
problem of the theory of electroelasticity, data on the time required to perform basic operations with stiffness matrices in 
various storage formats were summarized and presented in the form of a table. It has been established that the renumbering 
of grid nodes gives a significant increase in performance even without changing the internal structure of the matrix in 
memory. Taking into account the objectives of the study, the advantages and weaknesses of known matrix storage formats 
were named. Thus, CSR was optimal when multiplying a matrix by a vector, SKS was optimal when inverting a matrix. 
In problems with the number of unknowns of the order of 10°, iterative methods for solving a generalized eigenvalue 


problem won in speed. The performance of the software implementation of the Lanczos method was evaluated. The 
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contribution of all operations to the total solution time was measured. It has been found that the operation of solving 
SLAE takes up to 95% of the total time of the algorithm. When solving the SLAE by symmetric LQ method, the greatest 
computational costs were needed to multiply the matrix by a vector. To increase the performance of the algorithm, parallelization 
with shared memory was resorted to. When using eight threads, the performance gain increased by 40-50%. 

Discussion and Conclusion. The software modules obtained as part of the scientific work were implemented in the 
ACELAN-COMPOS package. Their performance for model problems with quasi-regular finite element grids was 
estimated. Taking into account the features of the structures of the stiffness and mass matrices obtained through solving 


the generalized eigenvalue problem for an electroelastic body, the preferred methods for their processing were determined. 


Keywords: piezoelectric materials, finite element method, sparse matrices, generalized eigenvalue problem, Lanczos 
method, Krylov subspace, preprocessing module, postprocessing module, triplets, coordinate storage format, compressed 


sparse row, CSR, compressed sparse column, CSC, Skyline 
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AnHoTayna 

Beedenue. LWinpoxkoe ucnomb30BaHHe Mbe3OMaTepHaJiOB B pa3IMYHbIX OTPacJIAX CTHMYJIMpyeT U3y4eHHe UX (PU3H4eCKHX 
XapaKTepHCTHK HU OOyCOBINBaeT AKTYaJIbHOCTb TaKHX W3bICKAaHHU. B paccMaTpHBaeMOM Cily4yae MOMAJIbHbIM aHasH3 
MO3BOJACT ONpeeIUTh padouylo yacTOTy HU KOIPPUUMEHT ZICKTPOMeXaHHYeCKON CBA3H MIbe3OIIEMCHTOB Pa3JIMYHBIX 
yCTpolcTB. ITH MHAMKATOPh! MpeACTaBIIALOT Cepbe3HbIi TeopeTHYeCKHH UM IpHkaqHOH uHTepec. Len uccneqoBaHua — 
pa3pa0oTKa YHCJICHHBIX MeTOAOB JIA peleHHA 3aa4H OlpeweeHuA YaCTOT pe30HaHca B CHcCTeMe yupyrux Ter. Ja 
OCTWKeHHA We HYKHBbI HOBbIC NOAXOAbI K WHCKpeTH3allMH 3aadH Ha OCHOBeE MeTOa KOHCUHBIX 3JIEMCHTOB 
BBINOJIHeHHe MpOrpaMMHOM peasH3alHu BbIOpaHHoro MeTOa Ha A3bIKe C# Ha IaTPopMe .net. AKTyaJIbHble pewieHHa 
CO3aHbI B KOHTeKcTe OuOMOTeKM KIaccoB KomiIekca ACELAN-COMPOS. OcnHopaHuple Ha oOpallleHHu MaTpHIl 
W3BECTHbI€ MCTOALI pewieHHA OOOOMIeHHOM 3aqa4H Ha COOCTBeHHbIe 3HaYeHHA HeEIPpHMeCHHMBbI K MaTpHuaM OosbUON 
pa3smepHoctu. Jia mpeoxomenua 3TOrTO OrpaHW4ueHHA B pecTaBueHHOH Hay4Ho padoTe peasM30BaHa sJIOrMKa 
MOCTpOCHHA MaTPHL Macc HM CO3aHbl IporpaMMHble MHTepdelich Wid OOMeHa TaHHbIMM O 3aa4ax Ha COOCTBeHHBIe 
3HadeHHA C MOJYIAMH Ipe- MW MOCTIpoLeccuura. 

Mamepuansi u memooo. J\na peanv3alyun 4HCIeHHbIX MeTOROB 3afelicTBOBaIH TaTpopMy .net HU A3bIK 
liporpamMMuposanusa C#. Barnyalua pe3syIbTaTOB HCCIeAOBaHUA NPOBOAMIACh IYTeM CpaBHeHHsA Hal eHHbIX 3Ha4deHHi 
C pelieHwiMu, ouyYeHHbIMH B W3BeCTHBIX CAE-nmakeTax (aHrs. computer-aided engineering — 
KOMIbIOTepH3HpOBaHHasd HHKeHepuA). Co3laHHble MOAMpOrpaMMbl OLCHMBAJIMCh C TOUKH 3peHHA MpOW3BOAMTeIBHOCTU 
MW MIpHMeHHMOCTH Ayia 3aqa4 OombuIoH pasMepHoctH. II[poBoqumMcbh 4YMCIIeHHbIe SKCIEPHMEHTH! C WesbIO Bas Waluu 
HOBBIX aJIFOPHTMOB B 3ajauax MasIOM pa3MePHOCTH, KOTOPbIe pellalOTCA W3BECTHBIMM MeTOAaMu B MATLAB. Jlanee 
TMOAXOA TecCTHpoBaIM Ha 3aqauax C OOJbIIMM UMCIIOM HeM3BECTHBIX HM C yYeTOM pactiapasWiesIMBaHHA OTCJIbHbIX 
onepauui. UroOn u30ex%KaTb Haxo7KeHHA OOpaTHOM MaTPHIbIl, IPOrpaMMHO peasM30BaIH MOAMPUUMpOBaHHbIM MeTOA 


Jlanyoma. Paccmotpemu dbopMaTbl xpaHeHHuA MaTpH B OMepaTHBHOM aMaATH: TpuIIeTE, CSR, CSC, SKyline. [na 
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pelieHHA CHCTeMbI JIMHeMHEIX areOpanyeckux ypaBHenuit (CJIAY) 3ayeticTBoBanM UTepallMOHHbI CHMMeTpHUHbI 
Metog LQ, afanTupoBaHHblii K ITHM (pOpMaTaM xpaHeHHaA. 

Pezyismamot ucciedoesanua. PaspaboTaHbl HOBbIe pacueTHbIe MOY, MHTerpupoBaHHble B OUOMMOTeKy KIaCcoB 
komiiexca ACELAN-COMPOS. IIpopegenbi pacueTsl AIA ompeyeeHHaA IpHMeHHMOCTH pa3IM4HBIX CopMaToB 
XpaHeHHA pa3spe2xKeHHBIX MaTPHI B ONepaTHBHOH MaMATH HU pa3sIM4HbIX MCTOAOB peaslu3alHu Onepalluit c paspexKeHHbIMU 
Matpuiamu. I padwyecku Bu3yaIM3HpoBaHa CTpyKTypa MaTpHL KECTKOCTH, MOCTPOeHHBIX [JIA OHO U TOM xe 3aa4H, 
HO C pa3yIM4HOH TepeHyMepalveli y3I0B KOHe4HOIeMeHTHOM ceTKH. IIpuMeHuTenbHO K 3aqa4ue TeopHu 
3JIEKTPOYIpyrocTH OOOOMIeHEI HU MpeACTaBIeHbI B Bue TAOHMIUbl TaHHbIe O BPeMeHH, HeEOOXOAMMOM Ha BBIMOJIHeHHe 
Oa30BbIX OMepalHii C MaTPHWaMH %KeCTKOCTH B pa3JIM4HbIX (POpMaTax xpaHeHHa. YcTaHOBIeHO, YTO TepeHyMepalHA 
y3J10B CETKH JaeT CyICCTBEHHbIM IpMpOcT MpOu3BOAMTeIbHOCTH JaxKe Oe3 H3MCHEHHA BHYTpeHHeli CTpyKTYpbl MaTpHUbI 
B TaMatTH. C yueTOM MOCTaBJICHHBIX 3aJad MCCIIeqOBaHHA Ha3BaHbl TIPeHMyIMecTBa UM cylaOble CTOPOHbI H3BECTHbIX 
cbopMaToB xpaHeHua MaTpuu. Tax, CSR onTuMasieH Ip yMHOXKeHHM MaTpHIbI Ha BeKTop, SKS — upu obpamenuu 
matpuusl. B 3aqjayax c YMCJIOM HEH3BECTHEIX TWopsKa 10° BEIMTpbIBaIOT B CKOPOCTH UTepaltMOHHble MeTOJbI pelleHus 
oOobmjeHHOH 3aqaun Ha COOCTBeHHbIe 3Ha4eHHA. OUCHHBaIaCh NPOM3BOAUTCIBHOCTb IporpaMMHOH peasM3alHu MeToa 
JIannoma. U3mepsiica Bkay BCex OMepalHii B OOMee BpeMa pelieHHa. BrlicHvsIOcb, YTO ONepalua pemseHua CJIAY 
3aHuMaeT 0 95 % oT oOmjero BpeMeHH paboTsI anropuTMma. IIpu pemenun CJIAY cummMetpwunbim MeToyom LQ 
HavMOOubUIMe BbINMCJIMTeIbHBIC 3aTpaTbl HYKHbI JIA YMHOXKeHHA MaTPHUbI Ha BekTop. Jaq yBemM4eHuA 
MpOH3BOAHTeIbHOCTH aIrOpHTMa UpHOermM K pacnapaswieMBaHHIo c OOmei MamaTBIO. [pH ucnoub30BaHHH BOCbMH 
MOTOKOB IIPOH3BOAUTEIbHOCTh BbIpoca Ha 40-50 %. 

O6cystcoenue u 3akmo4uenue. Ilonyuenuble B paMKax Hay4HOH paooTE! NpOrpaMMHble MOLYJIN OLIN BHEPeHbI B WaKkeT 
ACELAN-COMPOS. OneHeHa UX MTpOH3BOQMTeIbHOCTh JIA MOJCIbHBIX 3aad C KBa3HperyApHbIMU 
KOHCYHORJIEMCHTHBIMH CeTKaMH. C yyeTOM OcOOeHHOCTelM CTpyKTyp MaTPHIl 2%X€CTKOCTH HM Macc, MOyYaeMbIx pu 
pelleHuv OOOOMeHHOH 3aqaun Ha COOCTBEHHBIe 3HAYeHHA JIA 9IEKTPOYMpyroro Tea, OMpeeeHbI pe ANOUTUTeJIbHbIC 


MeTOJIBI JIA UX OOpadoTKH. 


KsroueBble CIOBa: Tbe30MaTepHalibl, MCTO, KOHCYHBIX IICMCHTOB, PaspexKeHHbIe MaTPHUbI, OOOOMIeHHaA 3ayayua Ha 
coOcTBeHHble 3HayeHHA, MeTog JlaHWoma, mosmpoctpanctBo Kppsiopa, MOfYIb lpempoleccuuHra, MOjyJIb 
MOCTIpOWeCCHHTa, TPHIWVIeTbl, KOOPAMHAaTHBIM (bopMaT xpaHeHHaA, CxKAaTHI pa3pexKeHHbIM pad, CSR, cxaroiii 


pa3pexeHHEI crombeu, CSC, SKyline 


B.saroqapHoctTH: aBropbl BbIpaxkaloT OarogapHoctp A.H. Conospesy u T.C. MapTsiHoBoli 3a NOMOLIb B pa3spaooTKe 
YMCJICHHBIX MeTOOB HU Poccuiickomy Hay4Homy PoHAy 3a PHHAHCOBYIO NOAepxKy UccIeqoBaHUA rpantom Ne 22—21— 
00318, https://rscf.ru/project/22-21-00318 


Aisa waTupopannsa. OraHecaH IT.A., UlreitH 0.0. Peanu3aima Oa30BbIxX oOlepallMii Wd pa3pexKeHHBbIX MaTpHIl B 
KOHTeKCTe pelieHHa OOoOMeHHON 3aqadH Ha COOCTBeHHbIe 3HadeHHA B KOMIIIeKce ACELAN-COMPOS. Advanced 
Engineering Research (Rostov-on-Don). 2023;23(2):121—129. https://doi.org/10.23947/2687-1653-2023-23-2- 121-129 


Introduction. Devices made of piezoelectric materials have been widely used, actively studied and improved for a 
long time. Medical ultrasound devices (diagnostic equipment, ultrasonic scalpels) [1-4] and mobile energy generators [5] 
should be noted separately. Paper [6] described the combination of photo- and piezoelectric effects to create efficient 
compact energy sources. New materials designed for the application under specific conditions are being studied in science 
and industry. In [7], the creation of a lead-free piezo-active composition suitable for operation at various temperatures 
was considered. 

In the study on piezoelectric elements, a key role is played by the modal analysis stage, which enables to establish 
the resonance and antiresonance frequencies of the device. These data: 

— are needed to find out the operating frequency of the device; 

— provide determining the electromechanical coupling factor — an important performance indicator of the device; 

— are input information in numerical experiments for problems on forced oscillations. 

The study was aimed at the creation of numerical methods for solving the problem of determining resonance 


frequencies for a system of elastic bodies. Achieving the stated goal required solving two tasks. The first was to develop 
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methods of discretization of the problem based on the finite element method (FEM). The second was to carry out a 
software implementation of the selected method in C# on the -net platform. All known programs take into account the 
context of the ACELAN-COMPOS class library [8]. When solving a generalized eigenvalue problem, methods based 
on matrix inversion are widely used. However, they are not applicable to large-dimensional matrices. In the presented 
research, this limitation was overcome as follows: 

— the logic of constructing mass matrices was additionally implemented; 

— software interfaces were created for exchanging data on eigenvalue tasks with pre- and postprocessing modules. 

Materials and Methods. Principally, the proposed approach was designed to solve static problems of electroelasticity 
when implementing the averaging method [9], which was used to calculate the effective properties of piezo composites. 
In this regard, only stiffness matrices were presented at the stage of constructing global FEM matrices. In this study, we 
additionally implemented the logic of constructing mass matrices and developed software interfaces (application 
programming interface, API) for exchanging data on eigenvalue tasks with pre- and postprocessing modules. The 
developed routines were evaluated in terms of performance and applicability for large-scale tasks. Numerical experiments 
were carried out to validate the algorithms created for such small-dimensional problems that provide obtaining a solution 
by general methods in the MATLAB computing package. Next, testing was performed on tasks with a large number of 
unknowns and taking into account the parallelization of individual operations. 


The mathematical model of the problem being solved consists of the defining relations [9]: 


PreO2U + pPpiau-V-o=$, V-D=0, (1) 

o=cf -(€+ Bs é)-ef -E, D+quD =e;- (€+Ga8) +95 -E, (2) 

6=(Vu+VuT)/2,E=—-Vo. (3) 
Here, o — stress tensor; p; — body density; ¢ — strain tensor; u — displacement vector; D — electric displacement 
vector; E — electric-field vector; f, — body force vector; » — electric potential; a,, By;, Sa — damping 


coefficients; cf, ef 


£. e™, 9% — tensors of elastic constants, piezoelectric modules and dielectric permittivity; index j — 


body number in the model. 
The discretization is performed by replacing: 
u(x,t)=Ni(x)-U(t), ot) = Nj (x): Od). 
Here, N,, — shape function matrix for the displacement field; Ny — shape function vector for electric potential; 


U (t), (rt) — global vectors of the corresponding nodal degrees of freedom. 


In this case, the original problem (1-3) takes the form: 


M-a@+K-a=F. (4) 
Here, matrices M and K are global matrices of mass and stiffness, respectively, and the vector is a general vector of unknowns: 
a=[U, ®]. 
In the problem of the theory of electroelasticity: 
0 K K 
M = M wa ). K _ uu ud . 5 
| 0 0 é uw —K ‘) ” 


Matrices M ws, Kw. and Ky — symmetric. In the case of harmonic oscillations at natural frequency @,, it is 
possible to write: 
a=v,sin(@,t), 
denoting the corresponding eigenvector by 1,. 
Consider free oscillations if F =0. In this case, task (4) is represented as: 
-—07 My, + K-v; =0. (6) 
Thus, the original problem is reduced to a generalized eigenvalue problem (6). For nonzero v,, inequality (6) is solved 


by finding the matrix inverse to K . However, at the same time, the sparse matrix becomes full, i.e., the method is 


unsuitable for large matrices. Therefore, it is needed to use other methods that do not require finding the inverse matrix. 
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To solve this problem, a modified Lanczos method was programmatically implemented in this paper [10]. The author of 
this modification is T. S. Martynova. The description of the development is not given in this article. Of the operations 
used in the method, the most expensive from the point of view of computational resources was the solution to a system 
of linear algebraic equations (SLAE), needed for performing a spectral transformation. 

Matrices M and K — sparse, with a small number of nonzero elements. Several formats are used to store such 
matrices in RAM: 

— triplets, or coordinate format; 

— CSR (compressed sparse row); 

— CSC (compressed sparse column); 

— Skyline storage format (SKS method). 

The coordinate format involves storing triples (triplets) of values (i, j, k), representing coordinates (i, j) and values (k) 
of nonzero elements. CSR is sometimes referred to as CRS or Yale format. It involves storing a sparse matrix in the form 
of three arrays. Consider matrix of N size with NZ nonzero elements. We describe the possible organization of its 
storage. All nonzero elements must be placed in one array of NZ size. The positions of these elements in the columns 
should be placed in another array of NZ size, and the third array of N size should be used to store the indices of the first 
elements of the rows. Similarly, the storage in CSV format is implemented. 

The SKS format assumes the storage of a variable-width matrix band that includes all nonzero elements. In this case, 
zeros are allowed. The efficiency of this format depends on the renumbering of the matrix rows. Methods for reducing 
the size of the tape are described in [11]; however, their applicability to the stiffness matrix obtained when solving a three- 
dimensional problem using FEM requires a separate study. 

To solve the SLAE (system of linear algebraic equations), an iterative symmetric LQ method (SYMLQ [12]) adapted 
to the storage formats listed above was used. 

Research Results. At the beginning of the study, we chose the optimal storage format for sparse matrices. The 
coordinate format enabled to quickly add and change an element of the matrix. These operations were needed at the stage 


of assembling the global matrix and taking into account the boundary conditions. In addition, for ill-conditioned matrices, 


to which K refers, a preliminary transformation is often used for normalization. It is also convenient to perform it in the 
coordinate format. However, this format is ineffective when it comes to algebraic operations. 

CSR is ill-adapted for changing the structure of the matrix: by adding a nonzero element, you need to insert into two 
arrays. In this case, the matrix is multiplied by a vector quite easily and efficiently. 

SKS has similar problems with the addition of nonzero elements and is highly dependent on the renumbering of 
unknowns in the problem. We focus on the example of a quasi-regular grid, which is used in the ACELAN-COMPOS 
package to work with representative volumes of composites. The width of the band containing all nonzero elements can 
be predetermined and depends on the number of nodes and the type of final element. In the general case of an arbitrary 
finite element grid, it is difficult to estimate the size of the band in advance. 

Four methods of numbering unknowns were used in numerical experiments. Figure | shows the structure of stiffness 


matrices constructed for the same task, but with different renumbering of nodes of a finite element grid. 
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Fig. 1. Stiffness matrix structure with various methods of node numbering: 
a — unknowns are ordered by nodes; b — first, displacement nodes are ordered, 
then potentials; c — nodes are sorted by layers of the FE grid, and unknowns — by nodes; 
b, d— nodes are sorted by layers of the FE grid, and unknowns — as in example 


Thus, the grid was a cube with regular partitioning by eight-node finite elements. A model matrix of 500 lines was 
used for the illustrations. The matrices shown in Figures | a and 1 b, were not subjected to additional renumbering of 
nodes and differed only in the numbering of degrees of freedom. In | a: 

a= [Hy My My Py yes Uyy sUyy Uy,» Py | < 

In 1 b, the unknowns responsible for the potential distribution were collected at the end of the vector: 

a= [Uy Myy Uy. yeep Uiyy My Uys > Pyro Py | : 

The unknowns in matrices 1 c and | d were numbered similarly, but the nodes of the finite element grid were pre- 
numbered according to their coordinates through alternately sorting all nodes by each of the coordinates. This technique 
is widely used to build more efficient SLAE solution modules, as it enables to work with the matrix in a suitable band 
format, convenient for parallelization. Similar external modules were implemented for the ACELAN-COMPOS 
complex [13-15]; however, in this work, only formats for storing sparse matrices of a general form were used. 


Table | summarizes the data on the time required to perform basic operations with matrices in various formats. 


Table | 
Time to perform basic operations with the stiffness matrix in the problem of the theory of electroelasticity. 
19,652 rows 
Storage . Elapsed time, ms 
Operation 
format la 1b lec ld 
CSR Conversion from coordinate format 123 132 97 117 
CSR Multiplication by vector, 100 operations 260 260 260 260 
SKS Conversion from coordinate format 690 703 124 268 
SKS Multiplication by vector, 100 operations 60,558 | 61,450 7,616 22,113 


The experimental results showed that the conversion operation from a coordinate storage format to a compact one 
took little time. At the same time, the renumbering of grid nodes to form a block-tape matrix made it possible to get a 
noticeable increase in performance even without changing the internal structure of the matrix in memory. CSR format 
turned out to be optimal in terms of the efficiency of the matrix-vector multiplication operation. When the matrix was 
inverted, SKS format was more efficient, but for problems with the number of unknowns of the order of 103, iterative 
methods for solving a generalized eigenvalue problem worked noticeably faster. 

Further, the performance of the software implementation of the Lanczos method was experimentally evaluated. The 
contribution of all operations to the total solution time was measured. As a result, it was found that the operation of 
solving SLAE took up to 95% of the total time of the algorithm. In the course of the algorithm, a Krylov subspace was 
constructed, and depending on its dimension, the number of SLAE that needed to be solved changed. Note that the 


Oganesyan PA, et al. Implementation of Basic Operations for Sparse Matrices when Solving a Generalized Eigenvalue 


dimension of the Krylov subspace was chosen based on heuristics with respect to the number of desired eigenvalues. 
Here, the SLAE differed only in the right part, so that the requirements for allocated memory remained low. Among the 
basic operations used in the course of solving SLAE by the SYMLQ method, the greatest computational costs were needed 
for multiplying the matrix by the vector. 

To increase the performance of the algorithm, the simplest parallelization with shared memory was implemented. 
Blocks of rows were allocated for the CRS format. They were transmitted to separate threads that calculated the 
corresponding components of the resulting vector. The performance gain was 40-50% when using 8 threads. At the same 
time, for matrices of the order of 10° elements, the increase was about 40%, and for matrices of the order of 10*— about 50%. 

Discussion and Conclusions. Within the framework of this study, a method for solving a generalized eigenvalue 
problem for matrices obtained by modeling electroelastic bodies was implemented. Software modules were created in C# 
for constructing mass matrices by the finite element method and performing auxiliary operations within the framework 
of the Lanczos method (working with Krylov subspace vectors, reorthogonalization, finding eigenvectors). The 
computational complexity was mainly due to the operations of multiplying sparse matrices by a vector. In this regard, 
numerical experiments were carried out to determine the optimal formats for storing matrices, the optimal structure of the 
matrix obtained as a result of renumbering the nodes of the FE grid and degrees of freedom in the nodes. A version of the 
SYMMLQ iterative algorithm using parallel computing was developed. The final scheme of work included three points. 
First, global matrices were constructed in coordinate format with a renumbering algorithm (Fig. 1 c). Secondly, the data 
was converted to CRS format. Thirdly, the data was processed by the Lanczos method, which included the SYMLQ 
method for solving SLAE. The results of the work were included in the ACELAN-COMPOS software package. 
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Abstract 

Introduction. The formation of the quality parameters of the surface layer and the operational properties of the parts 
occurs throughout all stages of their manufacture. However, the decisive impact is most often exerted by the stages of 
finishing. Therefore, in modern digital engineering, the task of process support of high quality of the surface layer of the 
part is one of the challenges in solving the problem of improving the quality and reliability and increasing the life cycle 
of manufactured machines. Surface plastic deformation treatment is instrumental in improving the performance 
characteristics of machine parts. Its essence is that the required quality parameters of parts are obtained not by removing 
a layer of material, but by plastic deformation. During the processing, both the dimensions of the parts and the physical 
and mechanical properties of the surface layers are changed. In this case, the technologist has the opportunity to 
significantly increase the life cycle of the manufactured products through controlling the process. These studies are aimed 
at providing the required quality parameters of the surface layer under processing with an eccentric hardener. 

Materials and Methods. The article presents the results of research on a new method of surface plastic deformation 
treatment — with an oscillating eccentric hardener. The considered processing method enables to obtain high quality of 
the treated surface, to process large-sized parts in places that are stress concentrators, to process welds, small areas of 
surfaces, whose hardening is needed for the part to fulfill its intended service. A set of theoretical studies was carried out; 
their results provided determining the parameters of a single interaction of the indenter and the surface of the part, the 
diameter of the plastic imprint and its depth. 

Results. Dependences for determining the surface roughness, the depth of the hardened layer and the degree of 
deformation were obtained. The resulting formulas were tested for adequacy by experimental studies. 

Discussion and Conclusion. The obtained research results can be used in the technological design of surface plastic 


deformation treatment processes. Further tasks for the study of the considered processing method are determined. 
Keywords: oscillating tool, eccentric hardener, surface roughness, hardened layer depth, degree of deformation 
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AnHoTanHa 

Beedenue. DopmMupoBaHue NapaMeTpoB KayecTBa MOBeEPXHOCTHOTO CJIOA HM 9SKCIVIyaTalMOHHbIX CBOHMCTB jeTael 
IIPOHCXOJUT Ha IPOT#KEHHH BCeX ITAMOB UX H3TOTOBJIeCHHA. OWHAaKO pelllaiolljee BIMAHHe Yallle BCETO OKa3bIBaIOT ITAIIbI 
(buHHmHOM oOpadotku. IlostToMy B cCOBpeMeHHOM LM@poBoM MallMHOCTpoeHHH 3ajaya TeXHOJOrM4ecKoro 
oOecreyeHHA BbICOKOTO Ka¥eCTBa MOBEPXHOCTHOTO CJIOA WeTaIM ABIIACTCA OHOM M3 BaKHeEMMMX MpH pelleHun 
TIpOOEMBI MOBBILNCHHA KAYeCTBA, Halex*KHOCTH HM YBeIMYeHHA 2%KM3HCHHOTO WHKa UPOM3BOAHMMBIX Mallu. Bexyuryro 
pOJb B HOBBIMCHHH 9KCIUIyaTalMOHHbIX XapaKTepHCTHK JleTalIeii MalIMH vrpaeT oOpadoTKa MOBepXHOCTHbIM 
TWIACTHYCCKUM esPOpMHupoBaHHeM, CYLWIHOCTb KOTOpOM 3aKIIOUaeTCA B TOM, YTO TpeOyeMble MapaMeTpbI KayecTBa 
WeTasleH JOCTHraloTca He yaIeHveM CJIOA MaTepHalia, a ero WlacTHYecKUM JedopmMupoBaHuem. B mporecce oOpaooTKu 
IIpOH3BOHTCA H3MeHeHHe Kak pa3Mepos JeTaliei, Tak MU (PH3HKO-MeXaHHYecKUX XapakKTepHCTHK MOBEPXHOCTHBIX CJIOeB, 
ylipaBlad KOTOPbIMH TeXHOJIOr HMCCT BO3MOXKHOCTL 3HAYHTeCIbHO YBCJIMYHBATb %KH3HCHHBIM WMKI MpOu3BOMMOl 
mpoayKuuuH. Llenbro HacToAmHx MccleqoBaHu sABIAeTCA OOecneueHHe HeOOXOJHMBIX TMapaMeTpoB KayecTBa 
MIOBCPXHOCTHOTO CJIOA Ip OOpadoTKe IKCLICHTPHKOBbIM YIIPOUHUTEJIEM. 

Mamepuansi u memooosi. B cratbe mpescTaBeHbl pe3ysIbTaTbI UCCIeqOBaHH HOBOrO MeTOa oOpadoTKH 
TIOBCPXHOCTHBIM IWIaCTHYeCKHM ePOpMHpoBaHHeM —  OCIIMJWIMPYOWJHM 9KCIICHTPHKOBbIM YIIPOUHHTesIeM. 
PaccmatpvBaeMbIli MeTO OOpaOoTKH MO3BONAeT MOyYaTb BbICOKOe KayecTBO OOpadboTaHHOM MOBepXxHOCTH, 
OCYIeCTBIATh OOpaboTKy KpylHoraOapHTHBIx JeTalielt B MeCTaX, ABJIAIOW[MXCA KOHI[CHTpaTopaMy HalipsKeHHH, 
oOpabaTbIBaTb CBapHBle WIBbI, HCOObIIMe y4acTKH MOBepXHOCTeH, yMpOuHeHHe KOTOPHIX HEOOXOAMMO JIA BbITIOHeHHA 
WeTasIbko cBoero cilyxKeOHOrO Ha3HayeHHA. BpinosHeH KOMIVIeCKC TeOpeTHYeCKHX MCCeqOBaHHH, 10 pesysIbTaTaM 
KOTOPBIX OMIpeylesIeHbl TapaMeTpbI eCAHHHYHOTO B3aMMOJeMCTBHA MHAeHTOpa C MOBePpXHOCTbIO JeTasIM, WMamerTp 
IWIaCTHYECKOFO OTTHeUaTKa HM ero rryOuHa. 

Pezyismamot ucciedoeanua. Wlonyuenbl 3aBMCHMOCTH JIA OpesesIeHHA WIepOxXOBaTOCTH MOBepXHOCTH, TtyOMHbI 
ympouneHHoro clot wv cTeneHu yedopmayun. Ilonyyenupre dopmysbl pou MmpoBepKky aeKBaTHOCTH 
9KCICPHMeHTAIbHBIMH HCCIeOBaHHAMH. 

O6cystcdenue u 3akmo4uenue. TlonyyeHHble pe3yiIbTaTbI MCCeqoBaHu MOryT ObITb UCHONb30BaHbI TIpv 
TeEXHOJOFHYeCKOM IpOeKTHPOBaHHU MpoleccoB OOpadoTKH MOBepXHOCTHBIM IWlacTHYecKUM JedopmMupoBaHHeoM. 


OnpezeseHsl WalibHenuine 3aadn TO HCCIICTOBaHHFO paCcCMaTpUBaeMOro MeTOa oOpadorTKn. 


KoroueBbie CJ10Ba: OocuWIMpyrounnt MHCTPpyMeHtT, 9KCLICHTPHKOBBIi YIpOUAHTeJIb, WepOXOBaTOCTb HOBepxXHOcTyH, 


rryOnHa yipOudHeHHOro CJIOA, CTCIICHb edopmaiuu 


BuaarogxapHocrn: aBTOPbI BbIPaxKaloT OlaroqapHOcTb pewakWHu *KypHaJia HW peweH3eHTaM 3a BHHMATCJIBHOe OTHOMICHHE 


K CTaTbe. 


Aisa waTnpopanna. Tamapkun M.A., Tamapxuu 9.9., Xamamt O. DopMupoBaHue Ka4ecTBa HOBEpXHOCTHOTO COM pu 
OTJeNOUHO-yMpouHsAoMeh OOpadoTKe MeTasei SKCICHTPHKOBbIM yupouHuTelem. Advanced Engineering Research 
(Rostov-on-Don). 2023;23(2):130—-139. https://doi.org/10.23947/2687-1653-2023-23-2- 130-139 


Introduction. The reliability and durability of machine parts largely depend on the quality of their surface layer. From 
numerous works on mechanical engineering technology, it is known that the formation of the quality parameters of the 
surface layer occurs at all stages of their manufacture. However, the dramatic impact is most often exerted by the stages 
of finishing. Therefore, in modern digital engineering, increased attention is paid to the production design of highly 
efficient finishing operations of parts, which enables to respond to a challenge of increasing their life cycle. Surface plastic 
deformation (SPD) treatment is critical in improving the performance of machine parts carried out under finishing 


operations. In contrast to traditional cutting methods, the quality parameters of the surface layer in SPD are obtained 
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through performing plastic deformation with special tools or working media. In the course of processing, simultaneously 
with the change in the size of the processed parts, the physical and mechanical properties of the surface layers are changed. 
In this case, the technologist has the opportunity to significantly increase the life cycle of the manufactured products 
through controlling the process. 

It should be noted that the widespread use of many SPD methods in industry is hindered by poor knowledge of their 
basic laws, difficulties arising in the process of designing optimal combinations of processing modes, and design 
parameters of the tooling. On numerous occasions, treatment modes are assigned on the assumption of the results of 
private experimental studies, which provides low processing efficiency [1—6]. 

The objective of these studies is to provide the required quality parameters of the surface layer during processing with 
an eccentric hardener. 

Materials and Methods. The need for applying SPD under the conditions of modern machine-building industries 
results in the creation of new processing methods. One of such methods is the processing of SPD with an oscillating 
tool — with an eccentric hardener. 

Figure 1 shows a kinematic diagram of an eccentric hardener consisting of vibrating body 1 suspended on flat 
springs 2. Vibrations of vibrating body 1, acting normally to the treated surface, are excited by the rotation of eccentric 
mass (unbalance) 3 around the vertical axis. The rotation axis of the eccentric mass is restricted from moving relative to 
vibrating housing |. The rotational motion is transmitted to the eccentric from electric motor 5 through flexible shaft 6. 
Tool head 4 with an instrument of the appropriate geometric form is attached to housing 1. The motion of tool 4 is limited 
by limiter 7 (a workpiece). In this case, the tool is an indenter with a spherical sharpening. It can be made in the form of 
a roller or a ball. The vibration system in the eccentric hardeners can be represented as a single-mass system with two 
degrees of freedom, which is under the action of a force varying according to the harmonic law. To study the system 
dynamics, we consider the features of its free oscillation under the action of centrifugal vibration excitation and the nature 
of the movement of the system hitting the limiter (a part). 

Under free oscillation, the vibrating system fixed at the end of flat springs 2 (Fig. 1) performs harmonic oscillations, 
which are excited by the rotation of eccentric 3 with a constant angular velocity. 

The proposed device can be effective when processing formed parts of not the most complex profile, and in some 
cases — when processing simple surfaces, such as planes or bodies of rotation. 

First, it is required to check the possibility of providing a wide range of the energy of the impact of the tool head on 
the surface of the workpiece in combination with relatively low altitude characteristics of surface roughness. 

Due to the lower rigidity of a flat spring (in our case, two springs) in the X direction, in comparison to the stiffness in 
the Y direction, the system describes a trajectory close to an ellipse with a larger semi-axis in the X direction. To analyze 
the law of motion of the system, we decompose the trajectories along the X and Y axes. 

The equation of motion of the center of gravity is nothing more than a mathematical expression of Newton's second law. 


Fig. 1. Scheme of eccentric hardener: 1 — housing; 2 — flat spring; 3 — eccentric mass; 4 — tool head; 
5 — electric motor; 6 — flexible shaft; 7 — limiter (workpiece) 
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Fig. 2. Rotations of the housing of the hardener during vibrations 


Research Results. To establish the main regularities of the process variables effect on the interaction of the oscillating 
indenter and the treated surface, it is necessary to take into account the kinetic energy of the indenter, the number of 
indenters, the radius of the indenter, the efficiency of the device, physical and mechanical properties of the 
workpiece [1-6]. 

Taking into account all the forces acting on the moving indenter, it is possible to write the equation of motion of the 


tool head in the Y direction: 


@ d 
m, y =-cy-p2+m,,,ro? COS Wf —™M, 8, (1) 
dt2 dt 
in the X direction: 
d2x dx : 
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where m, dt2’ mM, de projections of the inertia forces of the system on Y and X axes, respectively; ci, y, cx — 
ja t 
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projections of spring resistance forces on Y and X axes; uw ae uu = — projections of the resistance forces of the medium 
t t 


on Y and X axeS33) Mean @2 COS OF, Mean’ @2 Sin@t — projections of the perturbing force on Y and X axes; 


m.g — gravity (weight) of the vibrating system; m-— mass of the vibrating system; c1, c — spring stiffness in the Y and 
X direction; up — resistance of the medium; 77camn — mass of the eccentric; r— distance from the axis of rotation of the 
eccentric to its center of gravity; @ — angular velocity; t— current time value; y, x — current coordinate value. 

Due to the significantly greater stiffness of the springs in the Y direction relative to the stiffness in the X direction, 
the amplitude of the indenter movement in the Y direction is noticeably less than the amplitude in the X direction. 
Therefore, we assume that the system performs harmonic oscillations only in the X direction, i.e., we consider an indenter 
with only one degree of freedom. 

If we neglect the damped oscillation, then the equation of motion will have the form: 

x =bsin(at + £), (3) 


where b — amplitude of the oscillations; 8 — phase difference between the exciting force and the movements of the 
center of gravity of the indenter. 


Substituting this expression into equation (2), we find: 
2 


b = MNeam To 3 (4) 


y(c-o@'m,) +@°w? 


Machine Building and Machine Science 


133 


http://vestnik-donstu.ru 


134 


Advanced Engineering Research (Rostov-on-Don). 2023;23(2):130-139. eISSN 2687-1653 


Q, 
a (5) 
c-aom, 
Value u is determined from the expressions: 
o,m,d 
zee a 6 
On (6) 


where w; — frequency of natural oscillations; 5 — logarithmic decrement of attenuation. 
After differentiating equation (4) in time and examining the function at the extremum, we obtain an expression for 


determining the maximum speed of the indenter: 
v= Mag? O (7) 
V(c-o'm,) +@ Ww 


The largest kinetic energy of the indenter is determined from the equation: 


cam 


2 We-wmy+ow] 


2 2 2,6 
_mVy mm. O (8) 


The analysis of the interaction of spherical indenters and a deformable half-space (the surface layer of the workpiece) 
is described in the classical works of I. V. Kudryavtsev [1, 2, 4, 5]. The diameter of the plastic print can be determined 


from the dependence 


(9) 


In this case, the depth of the plastic print can be defined as: 


a ees (10) 
4\M-D,-HD 


where 7 — kinetic energy of the tool head; HD — dynamic hardness of the part material (ratio of the impact energy of 
the spherical indenter to the volume of the displaced material upon impact); D; — diameter of the indenter; n — efficiency 
of the device; M — number of indenters. 

When processing with an eccentric hardener, the roughness parameters of the treated surface can receive a constant 
(steady-state) value, which is reproduced during further processing of the surface of the part. The relief of the resulting 
surface can be both isotropic and anisotropic and is formed by repeatedly superimposing traces of a single interaction. 

When the oscillating indenter interacts with the initial protrusions of the microasperity, its height decreases with a 
simultaneous decrease in the depth of the cavities of the microasperity. With increasing processing time, the initial surface 
roughness profile is completely reshaped. As a result, a new microrelief is formed, and it has a specific character for each 
SPD method [7—19]. 

The finally formed roughness of the treated surface is called “established”. As a rule, its altitude parameters do not 
depend on the initial one. They are formed under the specific conditions of each processing method and depend on its 
process variables. Based on the methodology of papers [3, 4], a dependence was obtained for determining the steady-state 


surface roughness during treatment with an eccentric hardener: 


T: 
Ra =00075, |" (11) 
D,-M-HD 


The parameters of the surface layer hardening, which include the depth of the hardened layer and the degree of 
deformation, have a major impact on increasing the life cycle of the machined parts. As a result of theoretical studies, 


analytical dependences were obtained for their calculation during processing with an eccentric hardener: 


(12) 


(13) 
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The above dependences correspond to the physical meaning of the phenomena occurring under processing, and have 
been verified during complex experimental studies. 

In the course of the experimental studies, samples from various materials often utilized for the manufacture of machine 
parts were used: high-quality and alloy steels (steel 45, HVG, steel 30, steel 30HGSA, etc.), aluminum alloys (AL1, AVT, 
D16, etc.). Flat samples were treated with an eccentric hardener under different modes. Ball and roller indenters were used. 

According to the theoretical dependences, graphs of the dependences of the roughness of the treated surface, the depth 
of the hardened layer and the degree of deformation on the processing modes, characteristics of working media and 
processed materials were constructed. 

In the graphs (Fig. 3-6), a solid line shows curves constructed according to theoretical formulas, and the dots show 
the results of experimental studies. The construction of confidence intervals with a confidence factor of 95% has been 
performed. 

There is high convergence of the results, which indicates that the theoretical dependence reflects correctly the 


phenomena occurring during the processing of SPD with an oscillating tool — an eccentric hardener. 
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Fig. 3. Dependence of the surface roughness on the radius of the indenter: 
1 — steel 45 sample material; 2 — HVG sample material 
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Fig. 4. Dependence of the surface roughness on the Brinell hardness of the part 
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Fig. 5. Dependence of hardening parameters on Brinell hardness for various materials: a — on the depth of the hardened layer; 
b — on the degree of deformation 
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Fig. 6. Dependence of hardening parameters on diameter of the indenter: a — depth of the hardened layer; 
b — degree of deformation. 1 — steel 45 sample material; 2 — HVG steel sample material 

Discussion and Conclusion. Based on the results of the conducted research, the following conclusions can be drawn: 

1. A theoretical dependence has been obtained that provides determining the kinetic energy of the indenter under 
processing with an oscillating tool — an eccentric hardener. 

2. Dependences have been obtained for determining the diameter and depth of the plastic print, as well as the surface 
roughness according to parameter Ra, the depth of the hardened layer and the degree of deformation that provide 
predicting the quality of the treated surface. 

3. The results of theoretical and experimental studies of the treatment process with an eccentric hardener were 
compared. The discrepancy in the results did not exceed 15%. 

4. The process under study is subject to further investigation in order to determine other parameters of the quality of 
the treated surface, e.g., the magnitude of residual stresses in the surface layer, and [5, 6, 20] also to expand the range of 
treatment modes and design parameters in order to determine their optimal range. 

5. The dependences obtained for determining the key parameters of the surface layer quality make it possible to predict 


the results of processing and can be used to design processes with an eccentric hardener. 
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Abstract 

Introduction. In the areas of power engineering where the thermal energy of superheated steam is used, an important 
aspect of providing the reliability and safety of equipment is the heat resistance of the materials employed. In the 
manufacture of induction superheaters, the optimal material for the steam pipe (coil) is copper. However, its ultimate 
resistance to oxidation does not exceed 400 °C, which significantly limits the efficiency of steam generators. Therefore, 
the objective of the work was to study the kinetics of oxidation of the combined galvanic coating of the Mo-Ni-Cr system 
applied to copper tubular samples and intended for thermal protection of steam generator coils. 

Materials and Methods. A combined electroplating of the Mo-Ni-Cr system with a total thickness of 12-35 um was 
formed on the experimental copper tubular samples. A Mo sublayer with a thickness of about 1.5 um on the surface of 
the copper tube was formed to prevent the diffusion of Cu into the Ni coating. A 1.5 um thick chromium layer on the 
coating surface acted as an indicator of the oxidation process. A comparative analysis of the oxidation processes of the 
copper surface and the combined coating of the Mo-Ni-Cr system on a copper substrate was carried out using the methods 
of optical and electron microscopy, energy dispersive analysis, and precision determination of the growth parameters of 
oxide films. 

Results. The intervals of thermal stability of the copper substrate and nickel coating were experimentally determined. The 
obtained experimental dependences characterized the parabolic law of copper oxidation with the formation of a single- 
phase diffusion zone of CuO at temperatures above 350 °C, and nickel at temperatures above 750 °C, when the transition 
of NiO monoxide into oxide Ni2O3 began. The growth of oxide films according to quadratic laws provided a rapid increase 
in the thickness of the films, the accumulation of stresses in them, cracking, and chipping. 

Discussion and Conclusion. It is shown that the Mo-Ni-Cr electroplating is resistant to heating during long-term 
operation up to temperatures of 750-800 °C. The functional roles of Mo and Cr in the coating architecture were described. 
The work focused on the applied aspect of using the coating under study to increase the thermal stability of the steam 


pipelines of industrial induction superheaters with low and medium power. 
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AnHOoTanHa 

Beedenue. B tex oOnacTax 9HepreTHYeCKOLO MalIMHOCTPOeHHA, re HCHOMb3YeTCA TeNMOBad JHEPrua MeperpeToro napa, 
Ba)KHBIM aCIIeKTOM OOecre4eHHA HaexKHOCTH MU Oe30MaCcHOCTH OOOpyAOBaHHA ABJIACTCA TETIOCTOMKOCTb HCIOJIB3YEMBIX 
Matepuasos. IIpu v3roTOBsIeHHH HHAYKUMOHHBIX NaponeperpeBatesel ONTHMAJIbHBIM MaTepHasiOM JIA NaponpoBoya 
(3MeeBHKa) ABJIAeTCA Meb. OnHako e€ MpeyesbHad CTOMKOCTb K OKCHMpoBaHHIO He Mpespiuaet 400 °C, auto 
CyUlecTBeHHO orpaHHunBaert a:pdbeKTHBHOCTE paooTs! MaporeHepatopos. Ilostomy wWesbro padoTHI OblI0 HCcIesOBaHHe 
KMHeETHKM OKMCJICHHA KOMOMHMpOBaHHOrO TayIbBaHM4ecKOrO MOKpbITHA CHcTemMBI Mo-Ni-Cr, HaHeceHHoro Ha MeJHbIe 
TpyOuatTble OOpa3iibl W IpeqHasHayeHHOrO JIA TeMMOZaLIMTbI 3MCCBUKOB aporenepaTopos. 

Mamepuanvi u memoooi. Ha ompitHpix MefHbIX TpyOuaTBIX OOpa3zlax ObwIO cPopMHpoBaHO KOMOHHMpoBaHHoe 
rayIbBaHHdeckoe NoKpbiTHe cuctemp!l Mo-Ni-Cr c oOurel TommunHoli 12-35 mxm. Hoycnoi Mo tommnyol oKon0 1,5 MKM 
Ha HOBeEpXHOCTH MeqHOH TpyOKH OBI cPopMupoBaH AIA UpesoTBpawenua ZHpdpy3un Cu B Ni-noKxppiTue. Co xpoma 
ToMMMHON 1,5 MKM Ha NOBeEPXHOCTH NMOKPBITHA BbINOHA POJIb HHAMKaTOpa Wpowecca OKUCIeHHA. CpaBHUTesIbHbIi 
aHaJIHM3 MIpOWeccoB OKHCIICHHA MOBEPXHOCTH Me HW KOMOHMHMpoBaHHOrO HOKpbITHA cucteMbI Mo-Ni-Cr Ha MeqHoH 
TIOJJIOXKKE BBIMOJIHEH C MCHOIb3OBAaHHEM MeCTOAUK ONTHYECKON HM 3IEKTPOHHOM MUKPOCKOMNMH, IHeprogquciepCHOHHOTO 
aHasiH3a, a TaKoKe NIPeWM3HOHHOLO ONpeeIeHHA WapaMeTPOB poOcTa OKCHJHBIX IJICHOK. 

Pe3yiomamoel ucciedosanua. OKCTIepUMeHTANIbHO OMpeseeHbl WHTepBasIbl TepMH4eCKOM yCTOM4MBOCTH MeHOL 
MOJWIOKKH HW HAKeseBOro NOKpBITHA. Ilomy4eHHble IKCIePHMECHTAIIbHbIe 3ABHCHMOCTH XapakTepH3yI0T NapaoosMyeckHi 
3aKOH OKHCJICHHA Mea Cc OOpa30BaHveM OFHOa3sHon AN@dPy3HoHHon 30HbI CuO mpu Temmepatypax Bpitue 350 °C u 
HUKeJIA Ip TeMilepatypax Bbiie 750 °C, korga HaduHaeTca Mepexoy MoHooKcHya NiO nu B oxcugy Ni2O3. Poct okcHAHBIX 
TVICHOK MO KBapaTHYHbIM 3aKOHaM IIPHBOAUT K ObICTPOMYy YBCHYCHHIO TOJIMHHbI MWICHOK, HaKOIJICHHIO B HUX 
Hallps.KeHH, pacTpeCKMBaHHI0 H CKaJIbIBAaHHO. 

O6cyscdenue u 3aKkiio4enue. YloKa3aHO, 4TO TalbBaHwueckoe noKppitTHe Mo-Ni-Cr ycroitunpo kK HarpeBy TIpu 
JUIMTeIbHOM SKCIIyaTauHu BIVIOTb WO Temuepatyp 750-800 °C. Onncansr dyHKuMoHarbHbIe pom Mo wu Cr B 
apxHTeKType NoKpbITHa. Pabota akKWeHTupoBaHa Ha IIpHKa{HOM aciieKTe HCMOJIb30BaHHA HCCIeAyeMOro MOKPBITHA JIA 
MOBbIMMeHHA TepMuyeckol ycTowuMBocTu 3MeeBHKa-laponpoBoya TIPOMBILIJICHHBIX MHLYKIMOHHBIX 


llapolieperpeBatesieH Maso HW cpeqHel MOLIHOCTH. 


Kioueeoie_ cio6a: TlaporeHepaTopsl, TEIIOCTOHKOCTS, OKHCJIMTCIbHbIM Tipouecc, TaIbBaHHyecKHe IMOKpbITH4A, 


MUKPOCTpyKTypa, IIEKTPOHHaA MUKPOCKONMA, TpPaBHMeTpHYeCKHH aHasIH3 


Brazodapnocmu: aBTOPbI BbIpaxKaloT OnaroyjapHocTb petakWHu *KypHasia UW peWeCH3eHTaM 3a BHHMAaTeJIBHOe OTHOMCHHE 


K CTaTbe. 


Aaa waruposanna. Bapasxa B.H., Kygpaxos O.B., [punjenxo B.W. AcrtexTa! Term103aMTbI MaLIMHOCTpOMTeJIbHOTO U 
3HepreTHyeckoro OOopyAOBaHHA: MIpHMeHeHHe CTOMKUX K OKCHAMPOBaHHIO KOMOMHMPOBaHHBIX MOKPbITHM Ha OCHOBe 
Hukesia. Advanced Engineering Research (Rostov-on-Don). 2023;23(2):140—-154. https://doi.org/10.23947/2687-1653- 
2023-23-2- 140-154 


Introduction. Electrochemical deposition of metals is widespread in industry, being the basis of electroplating. One 
of the features of the development of this branch of science is usually attributed to the fact that its development took place 
“almost exclusively empirically” [1], starting primarily from the needs of various industries. Despite the fact that at 


present, the theoretical foundations of electrochemistry have been worked out quite deeply [2—6], the applied aspect is 
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still of priority importance here and determines most of the scientific tasks being solved, mainly related to the special 
conditions for the use of electroplating coatings. Under such particular conditions, superheaters are operated, for example, 
whose performance is associated with a significant change in the composition and temperature of steam along the length 
of the steam coil [7]. Depending on the power of the steam generator, the temperature of the coil along its length can vary 
from 150 to 650 °C, and for super-critical steam parameters in high-capacity modern steam turbines — even higher [8— 
10]. This work studies the possibility of using electroplating coatings to protect the steam pipe of an induction superheater 
from oxidation at high temperatures. On the totality of physical and technological properties (electrical conductivity, 
thermal conductivity, ability to plastic deformation, machinability by cutting, etc.), copper is currently an indispensable 
structural material for the manufacture of coils for household and low-power industrial steam generators. This is the 
reason for the interest in heat-protective coatings. However, the oxidizing ability of copper is also high, and its oxidation 
resistance does not exceed 400 °C. Based on the operating conditions of the superheaters under consideration, this 
circumstance offers a challenge of using coatings whose oxidation resistance level is above the thermal barrier of 600 °C. 

Materials and Methods. Taking into account the complex configuration of the coil, the presence of a large length of 
curved surfaces and its considerable overall dimensions, electrochemical deposition was chosen as the most 
technologically advanced method of applying a heat-protective coating. 

To select the composition of such a coating, accurate data on its operating modes were required. For this purpose, a 
thermal imaging analysis of the thermal operating conditions of an experimental induction three-coil six-turn steam 
generator with a capacity of 10 kW was carried out (Fig. 1) [11]. The steam pipe was made of profiled copper tube 
©25x1.5 mm of technical copper M2 grade according to GOST 617-2006. 


Fig. 1. Experimental induction superheater with copper steam pipe: a — general view; 


b — the most heavily oxidized sections of the steam pipe (coil) at the steam outlet 


Quantitative thermal analysis of the operating conditions of the steam pipeline, including forced operating modes, was 
performed using a no-contact thermal imager of Fluke Ti401 PRO model (manufactured by Fluke Corp., USA) [12] with 
the main technical characteristics: 

— infrared spectral range: from 7.5 to 14 um (long-wavelength); 

— thermal sensitivity: < 0.075 °C at an object temperature of 30 °C (75 mK); 

— error: +2 °C (at low temperatures) or 2%; 

— degree of protection: according to GOST 14254-96 (IEC 60529): IP54. 

The heat capacity of water vapor (c, = 33.6 J/(mol-K) under normal conditions) is approximately twice lower than the 
heat capacity of water (cp = 75.3 J/mol-K), which changes significantly the conditions of heat removal in the coil and 
contributes to the intensification of oxidation of the surface of the steam pipe; therefore, substantial heterogeneity of the 
degree of oxidation is observed along its length (Fig. 1 b). The thermal analysis results (Fig. 2) showed that the maximum 
heating temperatures were fixed at the output half-turn of the coil (Fig. 1 b). The range of their values was 530-540 °C 
with an absolute maximum of 541.38 °C. The average temperature values of most superheated (oxidized) half-turns were 
at the level of 420-460 °C. 

Taking into account the results obtained, in further studies on heat-resistant coatings to protect against oxidation of 


the coil surface, it is needed to focus on the maximum temperature load of 600 °C. In this regard, it is reasonable to use 
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nickel-based coatings. Ni forms the basis of most modern heat-resistant superalloys used in thermal power 
engineering [13-15], and the technology of electroplating Ni is quite well developed [15-18]. 

When applying experimental electroplating coatings to samples of copper tubes of technical copper M2, standard 
deposition modes and compositions of electrolytes containing Ni and Cr recommended by GOST 9.305 and 9.306 were 
used. During the operation of the steam generator, the coating is practically not subjected to mechanical action; therefore, 
it should not demonstrate outstanding mechanical properties during performance. At the same time, when applied to a 
curved convex surface, internal tensile stresses are formed in the coating. As the number of working heat shifts increases, 
their level grows. In this regard, the thickness of the investigated coating on the steam line should not be too large. It was 
taken as the average of the recommended in the literature ranges of values for nickel electroplating coatings performing 


protective and decorative functions, and was an approximate level of 20 um. 
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Fig. 2. IR image of a general view of the thermal fields of a working superheater (the window of the Fluke Connect program, 


the location of the coil is similar to Fig. 1 a, all values are in degrees Celsius) 


To study the microstructure of coatings and the topography of their surface, dual beam (electron/ion) scanning electron 
microscope ZEISS CrossBeam 340 (SEM) was used, which provides using a focused ion beam (FIB) to etch and perform 
cross-sections (sections of a given configuration) of samples directly in the vacuum chamber of the microscope with high 
positioning accuracy. The elemental composition of the studied surfaces was monitored using energy dispersive X-ray 
detector (EDAX) of X-Max 50N (Oxford Instruments) model, built into an electron microscope. 

The study of the oxidation kinetics of coatings was carried out through measuring the mass index of high-temperature 
gas corrosion (weight gain of the sample as a result of heating and oxidation). To determine the degree of oxidation of 
copper samples, both coated and uncoated, all samples were weighed before and after experiments at different stages of 
interaction. Gravimetric studies were carried out on analytical scales of “VLR-20” brand with a weighing accuracy of 10> g. 

Results and Discussion 

1. Qualitative analysis of oxidation kinetics. Properties of the chemical interaction of nickel and oxygen are 
manifested in the fact that Ni forms two modifications of monoxide: a-NiO with a hexagonal lattice (below 252 °C) and 
B-NiO with a face-centered cubic lattice. The transition occurs under continuous heating in the range of 250-300 °C. It 
was experimentally established that when heated to 630 °C, a diffusion process ran through a thin film of NiO monoxide, 
above 640 °C, a chemical process of NiO formation was established, which, when heated above temperatures of 800 °C, 
could cause the formation of Ni2O3 oxide [19]. 

The applied aspect of using a heat-resistant nickel coating on a copper substrate was complicated by two 
circumstances: the unlimited solubility of the Cu-Ni system components (Fig. 3) and the possibility of the Kirkendall 
interaction effect [21, 22] at the “Ni-coating — Cu-substrate” boundary. These phenomena reduced the resistance of the 


coating to oxidation due to the dissolution of copper in the coating. 
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Fig. 3. Diagram of the phase states of the copper-nickel system [20] 


To exclude the negative impact of these circumstances on the oxidation resistance of the coating, a combined 
electroplating of the Mo-Ni-Cr system was formed on experimental copper tubular samples (puc. 4). Mo sublayer with a 
thickness of about 1.5 um on the surface of the copper tube (Fig. 4 b) was formed to prevent the diffusion of Cu into the 
Ni coating under long-term operation of the steam pipe due to the practically insoluble system of Cu-Mo components and 
the limited solubility of Ni-Mo [20]. A layer of chromium 1.5 pm thick on the coating surface (Fig. 4 b) served as an 
indicator of the oxidation process (for more information, see below). The total thickness of the coating on the experimental 
samples with coatings was 12—35 um. The elemental distribution in the cross-section of the initial Mo-Ni-Cr coating 


(before the experiment with sample heating) is shown in Fig. 5. 


a) b) 
Fig. 4. Initial MoNiC coating in cross section, SEM: a — coating with thickness markers; b — homogeneous microstructure 
of the coating and hydrogen porosity at the boundary with the substrate 


Multilayer EMF 3 map Cr K series Cu K series 


3 

url 

3 

7) 

5 

no) = 2 ——— 

v 25um 25um 25um 
'S 

é a) b) c) 
3 

oO. 

P=) 

a 


144 


Varavka VN, et al. Aspects of Thermal Protection of Machine-Building and Power Equipment 
_—— 


Mo K series Ni K series 


ot —— 
25um 25um 


d) e) 
Fig. 5. Color maps of the distribution of chemical elements by the depth of Mo-Ni-Cr coating, EDAX: a — general view of the 


coating in cross-section (SEM); b—e — distribution of elements in general view image: Cr (b), Cu (c), Mo (d), Ni (e) 


For the experimental study of the oxidation kinetics under conditions as close as possible to the operating conditions 
of the steam pipe, the samples for the study were made of copper tubes with the corresponding wall thickness and 
diameter (Fig. 6). 

After plating with the technology that included heat treatment elements [11], the coating acquired a greenish hue 
characteristic of nickel monoxide NiO. When operating the steam pipe, the coating was applied only to the outer (convex) 
surface of the samples. However, in order to correctly determine the weight gain of the coating on the experimental 
samples, the coating was applied on both sides. 

Simultaneous heating of pure copper samples and coated samples was carried out at a fixed temperature in the range 
of 350-1000 °C in SNOL 6.7/1300 (2.4 kW) furnace in air. Exposure at a given temperature was 30 minutes. For the 
Statistical picture of the experiment, heating at each set temperature was carried out for S—7 samples with their separate 
loading into the furnace. The selective results of the experiment are visualized in Figure 7. 

From the experimental data obtained, it follows that under the conditions of the conducted heating, copper is relatively 
thermally-resistant to a temperature of 300-350 °C. At these temperatures, a dense thin oxide film of brown color is 


formed on the copper surface, regardless of its curvature (Fig. 7). 


Fig. 6. Samples for the study of oxidation kinetics; external surfaces of reference samples of pure copper (right) and copper coated 


with Mo-Ni-Cr (left), prepared for experiments 


e 750°C 
| | 650°C 


Fig. 7. Comparison of the outer surface of the samples after heating to the specified temperatures: on the left — coated 
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Fig. 8. Copper surface of the sample oxidized at a temperature of 650 °C, SEM: a — location of CuO crystals on the copper 


surface; b — morphology of CuO crystallites 


According to the literature data [23, 24], it appears to be two-layered: a thin sublayer of Cu2O is located on the surface 
of the sample, and a layer of CuO is located outside. Due to the small thickness of the latter, the internal stresses in the 
film are small. It has good adhesion to the substrate, low roughness, it does not loosen, and does not chip off the surface. 
As the heating temperature increases, the growth of the outer oxide layer CuO accelerates. Already at a temperature of 
450 °C, it becomes very loose and crumbles from the surface (Fig. | b). At the same time, a Cu2O sublayer of a 
characteristic reddish hue is found under it (Fig. 7). During further heating, peeling and shedding of the CuO oxide layer 
progresses — it practically does not stay on the surface up to temperatures of 650—700 °C. This, apparently, is due to the 
nature of crystallization of copper monoxide: its crystallites have a strict cut close to cubic (Fig. 8 b), weak conjugation 
with each other, and, most importantly, high heterogeneity of the nucleation sites (Fig. 8 a). Starting from the heating 
temperatures of 750 °C, the copper oxide film is compacted, the strength of its adhesion to the surface increases. On the 
concave surface of the samples, due to compressive configuration stresses, the CuO film is strong enough or can peel off 
completely from the entire surface of the sample without crumbling. At temperatures of 800-900 °C, the CuO oxide film 
behaves similarly on the outer (convex) surface of the samples. 

The surface of samples coated with Mo-Ni-Cr practically does not change up to a heating temperature of 750 °C. 
Upon further heating above 800 °C, the coating is first covered with a film of Cr2O3 oxide having a characteristic bright 
green color (Fig. 7). Then, at heating temperatures above 900 °C, the deeper layers of the coating are oxidized. Amphoteric 
chromium oxide Cr203 (Fig. 9 b) is fundamentally different in its morphology and crystallization character from CuO 
copper monoxide (Fig. 8). Cr2O3 oxide crystals have a characteristic polyhedron shape with a predominance of prismatic 
crystallites. Due to the large dispersion of crystallites in size, they, unlike copper monoxide, form a layer on the surface 


with a high packing density of crystallites. 
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a) b) 
Fig. 9. Surface of samples coated with Mo-Ni-Cr, SEM: 
a — after heating to a temperature of 650 °C; b — after heating to a temperature of 850 °C 
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If, prior to the oxidation of the chromium layer (below 800 °C), the coating surface has a very weak crystallinity 
character (Fig. 9a), then the appearance of Cr203 oxide gives the surface a well-known polycrystalline 
appearance (Fig. 9 b). In general, the nickel coating does not change its composition and structure by cross-section up to 
a temperature of 850 °C (Fig. 10). The chemical composition of the coating surface at this temperature indicates the initial 
degree of oxidation of the chromium layer. The presence of a thin layer of chromium oxide on the surface is confirmed 
by the data of energy dispersion analysis (EDAX) both by the depth of the coating and by the surface. Figure 11 shows 
that oxygen is concentrated in a much narrower surface layer (~1 pm) than the chromium layer (~3 um), which 


characterizes the initial stage of oxidation. 
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Fig. 10. Thickness and structure of Mo-Ni-Cr coating in cross section after heating up to 850 °C, SEM 


To reduce the volume of the article, the EDAX data is not provided in full. Its results show that the amount of 
oxygen on the surface increases from 30 to 50 at. % due to a similar decrease in chromium concentration, which 
indicates the oxidation of chromium (since the concentration of Ni does not change when heated from 650 °C to 
850 °C, and the presence of Ni in the detection results is caused by the penetration of X-ray radiation through a thin 
layer of chromium into the nickel base of the coating during EDAX analysis). The composition of the oxides 


corresponds to the compound Cr203. 


a) 


O, weight, % 
20 4 


16 = 
24 
8 
< 
0 10 20 30 40 50 hum 


b) 


Machine Building and Machine Science 


147 


http://vestnik-donstu.ru 


148 


Advanced Engineering Research (Rostov-on-Don). 2023;23(2):140-154. eISSN 2687-1653 
| 


Cr, weight,% 
30 = 
20 = 


10 2 


d) 


Fig. 11. Distribution of the main chemical elements in Mo-Ni-Cr coating by depth h after heating to 850 °C, EDAX: 
a — coating in cross section, the scanning direction (SEM) is shown; b—d — content of elements in the scanning direction; 


b — oxygen; c — chromium; d — nickel 


All stages of the oxidation process of the coating are shown in Figure 12. At temperatures above 800 °C, the surface 


of the coating starts to oxidize, as indicated by a change in its color — the coating acquires the Kensington Green color. 


a) 


Fig. 12. Successive stages of oxidation of nickel coating, optical microscopy, x100: a — initial stage: germination of nickel oxide 


(dark-gray precipitates) on the surface of chromium oxide (green field) at 850 °C; b — final stage: 
collective phase pattern in the area of coating chip at 1,000 °C, where: 


1 — chromium oxide Cr203; 2 — nickel oxide Ni2O3; 3 — copper surface (coating chip); 4 — sample edge 


This color corresponds to chromium oxide Cr203. Already at a temperature of 850 °C, rare single formations of Ni203 
nickel oxide can be found on the surface of the coating (Fig. 12 a, gray color), which, under further heating, gradually 
increase their area occupied on the surface. As they grow, which mainly spreads tangentially to the surface, under the impact 
of internal stresses, the Ni2O3 oxide film cracks and then chips off, exposing the surface of the substrate — pure copper of a 
reddish hue (Fig. 12 b). 

Thus, the basic result of this part of the research should be considered experimentally established temperature ranges 


of permissible use of materials for the manufacture of steam pipe of steam generator. Thus, a copper steam pipe without 


Varavka VN, et al. Aspects of Thermal Protection of Machine-Building and Power Equipment 


coatings is operable up to a temperature of 300 °C and can only be used to generate wet (not overheated) steam. Starting 
from a temperature of ~400 °C and up to ~700 °C, the CuO oxide film formed is very loose and easily crumbles from the 
copper surface. At higher temperatures, the film becomes denser, thicker, and its adhesion to the copper substrate 
increases. However, it is still prone to chipping during heat exchange. Its presence on the surface of the steam pipe 
significantly slows down the heat removal, and the chemical reactions of high-temperature gas corrosion that continue 
during heating, work in the direction of reducing the thickness of the pipe. Due to the heterogeneity of the ongoing 
processes, the operation of the steam pipe under these conditions becomes unpredictable from the point of view of 
emergency situations. The use of a combined electroplating of the Mo-Ni-Cr system increases the efficiency of the 
steam pipe to a temperature of 750-800 °C. When the temperature reaches 850 °C, the coating starts to oxidize along 
with copper. At 950 °C and above, the oxidized coating is prone to chipping, and its operation is subject to the same 
risks as the copper pipe. A distinctive feature of the investigated coating is self-testing: if the heating temperature 
exceeds 800 °C during operation, the surface layer of chromium turns the coating bright green and signals the danger 
of overheating. The indicator layer of galvanic chrome after oxidation can be easily restored, and the operation of the 
steam pipe then continues. 

Quantitative analysis of oxidation kinetics 

As part of the performed studies, a quantitative analysis of the kinetics of oxidation of pure copper samples and 


Mo-Ni-Cr coated samples was carried out. The specific mass gain MM = Am/S , observed during the heating process, 


was used as a measured parameter, where Am — increase in the mass of the sample, g; S — area of the oxidized surface 
of the sample, cm”. According to the qualitative analysis method described above, the experimental data of M values 
were obtained during the heating of tubular samples made of pure copper and coated samples. They are presented in 
Tables | and 2. 

The Tables show the spread intervals of the data obtained for fixed values of heating temperatures (Table 1) or the 
holding time in the furnace (Table 2), as well as the average value of M from each interval. 

Statistical processing of the data shown in Tables 1 and 2, performed using the MathCAD application software 
package, which included interpolation procedures, allowed us to obtain kinetic dependences shown in Figures 13 and 14. 
The rectilinear graphs of the obtained dependences, shown in Figure 13 b in Arrhenius coordinates (—In M — 1,000/T), 
characterized the parabolic law of oxidation of copper at temperatures above 350 °C and nickel — at temperatures 
above 750 °C [19, 25-27]. 


Table | 
Experimental data on the specific mass gain of samples under furnace heating in air for 30 min 
No. of Furnace Specific weight gain M, 10° g/cm? 

experiment temperature, °C Uncoated copper tube Ni-coated copper tube * 

1 350 1.35 + 0.31 = 

2 450 6.23 + 1.37 = 

3 550 24.73 + 3.60 a 

4 650 58.51 + 9.02 1.52 + 0.24 

) 750 138.88 + 17.85 4.25 + 0.58 

6 850 241.03 + 25.25 13.12 + 1.08 

i) 1,000 564.70 + 49.76 35.71 £2.78 


* minimum weight gain equal to 10> g, measured on the analytical scales, was taken for the criterion of the absence 


of oxidation (dash in the Table) 
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Table 2 
Experimental data on specific mass gain of samples at different exposure time in the air of the furnace 
No.of Holding time in the Specific weight gain M, 10> g/cm? 
experiment furnace, min Uncoated copper tube at 600 °C Ni-coated copper tube at 800 °C 
1 5 13.68 + 2.02 3.03 £0.41 
2 15 29.43 + 4.15 5.93 + 0.54 
3 30 35.76 + 4.84 7.11 £0.67 
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Fig. 13. Temperature dependences of mass gain M of pure copper samples (1) and samples coated with Mo-Ni-Cr (2): 
a — in absolute units; b — in relative coordinate system 
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Fig. 14. Kinetics of time variation of mass gain M of pure copper samples at 600 °C (1) and samples 
coated with Mo-Ni-Cr at 800 °C (2) 
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The growth of oxide films according to quadratic laws occurs with the formation of single-phase diffused zones, in 
this case consisting of CuO and NiO oxides, respectively. It causes a rapid increase in the thickness of the films, the 
accumulation of stresses in them, cracking and chipping. An additional contribution to the acceleration of this process 
is made by the curved outer surface of the copper tube [28, 29]. 

Conclusions 

1. The performed set of studies has shown that the combined electroplating of the Mo-Ni-Cr system is a 
sufficiently effective protection of the copper steam pipe from oxidation. The coating is able to provide a long-term 
operation of the steam generator up to heating temperatures of 750-800 °C. 

2. Long-term heat resistance of the coating is provided by an internal Ni layer with a recommended thickness of 
20-30 um. The study of the oxidation kinetics of the coating, performed by optical and electron microscopy, energy 
dispersion analysis, as well as using precision methods for determining the growth parameters of oxide films, has 
shown that the nickel coating is indifferent to heating up to temperatures of 600—650 °C. In the temperature range 
700-900 °C, the oxidation of the coating occurs with the formation of NiO monoxide according to the parabolic law. 
At higher temperatures, oxidation progresses due to the formation of Ni2O3, oxide film, which quickly causes its 
growth, cracking and chipping. 

3. The combined architecture of the investigated nickel coating includes two thin layers of Mo and Cr. The Mo 
sublayer with a thickness of about 1.5 um is located on the surface of the copper tube (substrate). Its function is to 
prevent the mutual diffusion of Ni and Cu at the coating — substrate interface during long-term operation of the steam 
generator, since the dissolution of copper reduces the heat resistance of nickel and disrupts the performance of the 
coating. The outer layer of chromium with a thickness of 2—3 um serves as an indicator of the degree of oxidation of 
the coating. The first sign of excessive oxidation of the coating is the appearance of a bright green hue on the surface 
of the coating, which is associated with the formation of chromium oxide Cr2O03 at temperatures >800 °C. The 
overheating indicator — a layer of chromium — is easily updated and contributes to the prolongation of the life cycle 


of the steam generator. 
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Abstract 

Introduction. The assessment of the manufacturability of products — as a stage of production planning and a key aspect 
of the development of modern industrial machining systems — is an urgent task of modern mechanical engineering. In 
this regard, theoretical and practical research on the development of methodological approaches to determining the weight 
significance of quantitative indicators in assessing the manufacturability of parts is highly relevant. The objective of the 
presented work was to develop an evaluation method aimed at improving the quality of part processing and the 
effectiveness of the performance of multiproduct manufacturing systems based on the development of additional 
quantitative indicators for assessing production manufacturability. 

Materials and Methods. To assess the impact of quantitative production indicators associated with time spent during 
equipment downtime, a model was created. It was aimed at predicting event flows of delivery of batches of parts for 
manufacturing for a specific operation and flows of processed parts using the queuing theory apparatus. This approach 
makes it possible to take into account both the design-engineering characteristics of parts, the features of a particular 
production system, and the emerging manufacturing situation. 

Results. The degree of influence of the manufacturability indicators at the level of the process operation was determined 
by assessing the possible impact on the components when calculating piece-calculation time (Tim.x.). The interrelations 
between the manufacturability indicators and expenses for all items of the production cost of part processing (Com), as 
well as costs associated with organizational downtime of equipment (Cyy...i) were established. The degree of influence of 
the indicators of manufacturability relative to other indicators was determined by using the apparatus of paired 
comparisons in decision-making in relation to all structural elements of production costs. 

Discussion and Conclusion. The approach to the implementation of this design procedure was described, which provided 
taking into account the composition and capabilities of processing equipment of a particular production and the actual 
production situation. The developed formalized models make it possible to comprehensively predict the impact of the 


manufacturability indicators of parts on the performance effectiveness of machining systems during their manufacture. 


Keywords: production planning, product manufacturability assessment, quantitative indicators of manufacturability, 
machining production systems, production efficiency 
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AHHOTalna 

Beedenue. Owenka pov3BOACTBeHHOM TeXHOJIOTMYHOCTH V3rOTaBJIMBaeMBIX W3qeuMM — 9Tall TexXHOJOrM4ecKoH 
MOATOTOBKH HM KJIIOUeBOM aCIIeKT pa3BHTHA COBPCMCHHbIX MPOH3BOACTBCHHBIX MexaHooOpabaTHIBaIOLIHX CHCTeM — 
ABJIACTCA AKTYAJIbHOM 3aaueli COBPeMeHHOTO MallMHOcTpoeHus. B sto cBa3H TeopeTH4ecKHe MH TpakTH4ecKHe 
MCCIeqOBaHHA TO pa3paOoTKe MeTOAMYeCCKHX MOAXOAOB K OMpefeseHHIO BeECOBOM 3HAaYMMOCTH KOJIMYCCTBCHHBIX 
ToKa3aTeleH Mp OWeCHKe MpOH3BOACTBeEHHOM TeXHOJIOrHYHOCTH JeTaIeli ABIIAIOTCA BECbMa aKTYaJIbHbIMH. Llemb10 
TIpesCTaBIeHHOM paOoTEI ABMJIACbh pa3paOoTKa MeTOJ{a OL[eCHKH, HaleIeHHOro Ha MOBbIMIeHHe KayecTBa OOpadoTKH 
WetaseH wu 3xpdexTuBHocTH PyHKUMOHHPOBaHHA MHOTOHOMEHKIaTYPHbIX TIPOH3BOACTBCHHBIX CHCTeCM Ha OCHOBe 
pa3paOoTKH JONOJHUTeIbHBIX KOJIMYCCTBCHHBIX HOKa3zaTeel OMCHKM MpOH3BOACTBCHHOM TeXHOJIOrM4HOCTH. 
Mamepuaavi u memoovi. [na oueHkKu BIMAHHA KOJIMYCCTBCHHBIX MPOM3BO]CTBCHHBIX OKa3aTeslel, CBA3AHHBIX C 
3aTpaTaMH BpeMeHH Ip Mpoctoe obopyAOBaHHA, CO34aHa MOJ{eJIb IPOTHO3HPOBaHHA MOTOKOB COOBITHH MocTyMIeHHA 
TapTHi etTasei Ha M3rOTOBJIeHHe Ha OMpeeIeHHY!O OepallHio HW MOTOKOB OOpaoboTaHHBIX JeTalei C UCHOsb30BaHHeM 
alllapatTa TeOpHi MaccoBoro OOcIYKUBAaHUA. Tako MOAXO MO3BOIACT YYCCTL, KAK KOHCTPYKTOPCKO-TexHOJIOrHyecKHe 
XapaKTepHCcTHKH jleTalieii, OCOOCHHOCTH KOHKpeTHOM MpOM3BOCTBeHHOM CHCTeMbl, TaK HU CKJIa{bIBalOUyIoca 
IIpOH3BO]{CTBeHHYIO CHTyallHto. 

Pezyismamot uccaedosanua. WocpexCTBOM OW€HKH BO3MO2KHOFO BIIMAHHA Ha COCTABIAIOWIMe IIpH pacueTe WITydHO- 
KaJIbKYIAWMOHHOTO BpeMeHH (Tiwmx.) Ha YPOBHe TeXHOJOrM4eCcKON onepallun Obvia oMmpeyeweHa CTeHeHb BIIMAHHA 
Toka3aTejieH TEXHOJIOTHYHOCTH. YcTaHOBJICHbI B3AHMOCBA3H Me*K Ay WOKa3aTeIAMM TEXHOJIOPHYHOCTH HM 3aTpaTaMH 110 
BC€M CTaTbAM TeXHOJIOrM4ecKOH cebecTOHMOCTH OOpaboTKH 3aroTOBKH (Coy), a Take 3aTpaTaMU, CBA3AHHBIMH C 
OpraHv3al{HOHHBIMH TIpoctosamMu OOopyAoBaHHA (Crpo.i). C MOMOLIbIO IIPHMeHeHHA allliapatTa MapHbIxX cpaBHeHuii pu 
IIPHHATHH pelWieHHH MpHMeHHTebHO KO BC€M CTPYKTYPHbIM 9JIEMCHTaM IIPOM3BOACTBeCHHBIX 3aTpaT onpeyemena 
CTeIHCHb BIMAHHA MOKazaTeei TEXHOJIOFHYHOCTH OTHOCHTEIbHO APyrHx WoKa3zarTesen. 

O@6cystcdenue u 3aKiio“enUue. OrmcaH MOAXOA K BBIMOJHeEHHIO WaHHOM MpoeKTHOM Mpouelypbl, NO3BONAIOWIMH 
YUYHTEIBATh COCTaB H BO3MO2%KHOCTH TeXHOOrMYecKOrO OOOpyAOBaHHA KOHKpeTHOrO MpOv3BOACTBa HM pealIbHO 
CKabIBaIOWyIoca IpOM3BOACTBeHHy!o cuTyauio. PaspadoTaHHbie dopMayIH30BaHHbIe MOJeIM MO3BOIAIOT 
KOMIVICKCHO CIIPOrHO3HpoBaTb BUIMAHHe ©oKa3aTeelH + TeXHOJIOrHYHOCTH eTalelH Ha 93Pd@eKTHBHOCTB 


(byHKUMOHHpoBaHHA MexaHooOpaOaTbIBarlollHx CHCTeM Hpi UX H3PrOTOBJICHHH. 


KoroueBbie CI0Ba: TeEXHOOrM4ueckaar MOATOTOBKa IIpOH3BOACTBa, OWlCHKa TCXHOJIOTHAYHOCTH v3 qe0Mi, KOJIUUCCTBCHHBIC 
TlOKa3aTeJIA IIpOu3BO]CTBeHHON TeXHOJIOTHYHOCTH, MexaHooOpaOaTBIBarolle TIPOH3BOACTBCHHbI€ CHCTCMBI, 


addekTuBHocts (PyHKUMOHHpOBaHHA IPOM3BOCTBAa 


Baarogxapnocrn: aBTOPbI BbIPaxKaroT OlaroqapHOcTb peqakWHH U peleH3enTy 3a BHHMAaTeJIBHOe OTHOMCHHE K CTATbe HU 


BbICKa3aHHble MIpeAIOKCHHNA, KOTOPbIe MO3SBOJIMJIN MOBbICHTb Ce KadyeCTBO. 


Aaa waruposanua. bouxapes IL.1O. Kopones P.J]., boxospa JI... Komnmexcuaa oeHka mpov3B0CTBeHHOM 
TexHOOrM4HocTH u3qenul. Advanced Engineering Research (Rostov-on-Don). 2023;23(2):155—-168. 
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Bochkaryov PYu, et al. Comprehensive Assessment of the Manufacturability of Products 


Introduction. The development of machine-building production under modern conditions is impossible without a 
serious increase in scientific research related to the development of theory and methodological principles of formalization 
of all stages of product production, which are the basis of future intelligent support systems for the creation and 
manufacture of technical objects. In this aspect, the solution to the tasks of design and process planning is a challenge [1, 
2]. Despite numerous works in this area, automated systems that provide for even minor functional actions of designers 
and technologists related to the implementation of creative design solutions have not been created yet. 

A prerequisite for the production planning of the effective functioning of machining systems is monitoring and 
analysis of the current production situation, as well as information about the condition of equipment and engineering 
support. Rational production decisions can be made only based on full knowledge of the above. Even an experienced 
technologist is not able to collect and analyze such a large amount of information. Therefore, decisions are often made 
subjectively and unreasonably, the design of processes and their implementation are spaced out in time, and the use of 
computing systems is hindered by the lack of models describing the process of production planning. 

R&D works on the creation of a system for planning multiproduct processes are devoted to solving the tasks 
formulated [3]. They are based on a conceptual approach to the formalization of all design procedures for providing the 
process planning of machining industries, taking into account specific features, capabilities of equipment and tooling. 
One of such design procedures is the assessment of the manufacturability of products, which is traditionally given 
insufficient attention. The role of this stage is significantly underestimated. 

All performance indicators of the production system operation are determined by the complexity of the products and 
the degree of production capacity. There is often inconsistency between these two indicators, which causes the inability 
to meet the requirements for the quality of products, downtime, and irrational use of equipment. Objective data on the 
feasibility of manufacturing products in a specific production system, along with known tasks and methods of solving 
them, should be obtained precisely when evaluating the manufacturability of products. 

Materials and Methods. Scientific studies on creating formalized models for establishing links between engineering 
and design tasks for the preparation of industrial machine-building systems are of great importance. Due to the increasing 
global rivalry in the manufacturing sector, the primary task is to increase the efficient operation of equipment during the 
implementation of production processes, taking into account compliance with the specified requirements for the quality 
of parts, which, in turn, are installed during the design process. 

The challenges of the modern conditions of the operation of industrial complexes involve providing the 
manufacturability of products. Currently, methods for assessing the manufacturability of products, taking into account 
the need for compliance with the requirements of standards, directly depend on the qualification of the technologist 
(designer) and their knowledge. This approach does not guarantee making reasonable engineering decisions and hinders 
the automation of project procedures. 

The assessment of manufacturability as a stage of pre-production is carried out to establish the relationship between 
the costs of manufacturing the product and its design features. The results of such an assessment are often contradictory, 
there is no complete mathematical description of the procedure for its implementation. 

To resolve the current situation, it seems appropriate to implement the following steps in practice [4-6]: 

— establishment of relative weight characteristics of manufacturability indicators based on the parameters of products. 
The solution to this problem at the stages of development of working design documentation, when there are no engineering 
solutions for manufacturing, is difficult to implement, but paramount; 

— development of the existing range of quantitative indicators for the implementation of the procedure for assessing 
manufacturability; they should provide taking into account specific approaches to the production planning for particular 
industrial complexes. 

The creation of methodological support for the design procedure of assessing the manufacturability of products should 
be based on an extensive design and engineering database that takes into account its structure and the relationship between 
the elements of models used in the design and implementation of the processes. The planning system of multiproduct 
processes meets these requirements and enables, along with the possibility of evaluating known and used quantitative 


indicators in production, to offer new ones [7]. 
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In accordance with the principles laid down when creating a process planning system, the key performance criterion 
is the operating time of the production system for manufacturing the whole set of products. It includes all the costs of the 
production cycle and is directly related to the cost of production of parts. Given this situation, the authors propose an 
approach that provides a conclusion about the significance of these indicators for specific production conditions. The 
approach is based on the establishment of relationships between the elements included in the estimation of the cost of 
manufacturing parts, and quantitative indicators of manufacturability. 

Research Results. The sequence of implementation of the developed approach includes several design procedures 
that take into account both the design features of the parts being processed and the organizational and technological 
features of the production system, including the composition and capabilities of the equipment, as well as the specifics of 
the program of manufactured products. 

Initially, at the level of the process operation, the degree of influence of the manufacturability indicators was 
established by assessing the possible impact on the components when making the time per piece calculation (Tim.x.). 
Figure | shows the structure detail (Tium.x.) for the turning operation, through which the analysis was carried out and the 
possibilities of the impact of production performance indicators on each individual value in the calculations (Tym...) were 
established. Similar studies, which provide establishing analytical dependences between quantitative indicators of production 


manufacturability and structural elements of process operations, were performed for other groups of process facilities. 


yem.cu + r. M.3.0 + L ses 4. pew. + Ls cca T: / er ate L seam pu. or: L sera tise. + rT, c.y. ot Li, npoe. a7 L sein: F T a6 npoe. + D vem.XZ a La m.COK 
‘i H oF ls 4 + L ina. + La + Lapa u + Li, Ku + l, np + / ee ay T, pao. + onep. + 0.mex.OOK a ocM 342 + uncmp.macm. 
ad 
T T i or ‘ia F Tots 
= nap |_ 
Fone =| ——— |= 2 aan t+ 
n,, N,, 
eA Ox f 
I+, 41 
1 2 3 
ia + rn + ifm ‘i XYZ + (ne 1M + T aieseuy 1M + Li dee a Lp pos: + scat ai he 
L + 
hie ix 
ji nxs <— T 47 T +T 
o 6 mex.00c ope 00€ 
1000xv S,, 
1x d n, ja + Tis F : ae 


| \ 


oor 
nep. pao yem.cu 
v,xK,xK,xK,xK,.xK; ihe ; 
u Y 1) gl 6 


Fig. 1. Block diagram Twm.x. under turning 


The interrelations between the manufacturability indicators and the expenses for all items of the production cost of 
machining workpiece Coy (Fig. 2), as well as the expenses associated with organizational downtime of equipment, are 
established Cnp.o.i.(Fig. 3). 
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To assess the impact of quantitative production indicators associated with equipment downtime, a model has been 
created for predicting event flows of delivery of batches of parts for manufacturing for a specific operation and flows of 
processed parts using the queuing theory apparatus. This technical approach is used in the process planning system [8, 9]. 
As an example, Figure 4 shows the results in the form of a Gantt chart. This approach enables to take into account the 
design-engineering characteristics of parts, the features of a specific production system and the emerging production 
situation. 

Figure 5 shows an enlarged diagram of the structure for determining production costs under manufacture of parts, 
used to assess the specific weight of quantitative indicators for assessing manufacturability in the process planning system. 
The analysis of the possibility of the influence of each indicator on the efficiency of the entire production system in the 
manufacture of a batch of selected parts for specific production conditions was carried out. 

The results of the presented analysis and the established relationships between the manufacturability indicators and 
the efficiency of machining systems allowed us to move on to solving the issue of establishing the significance of 
quantitative indicators of industrial manufacturability. The presented fragment (Fig. 6) contains information about the 
above links in relational form and is supplemented with information about the specific weight of cost elements (as a 
percentage). The data are obtained on the basis of statistical processing of the results of the real production system 


operation. In the absence of such information, it is possible to use general machine-building or industry-specific data. 
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Fig. 6. Fragment of the impact of quantitative indicators of production manufacturability on cost components 


To establish the weighting significance of the indicators of production manufacturability, it remains only to solve the 
problem of determining their impact directly on each element in the presented structure of production costs. The solution 
to this problem was carried out using the apparatus of paired comparisons in decision-making (the iterative Berge’s 
process [10]). This method provided determining the degree of impact of manufacturability indicators with other 
indicators in relation to all structural elements of production costs. Table 1 shows a pairwise comparison of technological 


performance indicators relative to the basic time as an example (To). 


Table 1 
Pairwise comparison of manufacturability indicators relative to basic time (To) 
To (14 %) 
Measured i) elle hae| 4 = 4s 
indicators o . 3 |}4)]5 ]/6]7 P| a se 23 2 2 : 
11 — 1);0;,0;0), 1 2 1 2 2 13 0.098485 
12 0 = 1);0;}0;2)],0) 0 1 1 2 2 9 0.068182 
13 0 1 -/|/0O]1}]1 707; 0 2 2 0 0 7 0.05303 
14 1 2 2/-]1)]14) 1 1 1 1 1 1 13 0.098485 
15 2 2 2)/1)-) 17) 1 1 1 1 2 2 16 0.121212 
16 2 2 1); 1}]0;-]0) 0 1 1 0 0 8 0.060606 
17 2 0 1; 1}] 17 ty-) 1 1 1 0 0 9 0.068182 
19 0 2 2/1} 17,2 )2 4) - 1 1 0 0 12 0.090909 
20 0 2 oO; 1}, 17141 1 = 1 1 1 10 0.075758 
21 1 1 Oo; 1)]1)17i1 1 1 = 1 1 10 0.075758 
22 0 0 2/1}),2;);2)2) 1 1 1 13 0.098485 
23 0 0 2/1 }]0;2 4,24) 2 1 1 1 — 12 0.090909 
Sum | 132 1 


Industrial testing and approbation of the developed models was carried out under the conditions of 
“GAZPROMMASH” LLC, specializing in the batch production of direct-acting gas heaters with an intermediate coolant 
and a modified series of stations, regulators, filter blocks and valves, high-pressure valves. The initial data for the 


experiments were: a generated and completed database containing information on the process capabilities of the 
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equipment (production unit No. 1), information on the actual condition and technical and economic characteristics of the 


site performance (Tables 2-6), the program of manufactured parts (drawings of individual parts are shown in Fig 7). 


Table 2 
Ratio of the components of production costs 
(site no. 1, LLC “GAZPROMMASH” LLC) 
306u. = 5,682 (100 %) 
Con = 4,662 Rub. (78 %) Cnp.o.i = 1,020 Rub. (22 %) 
= x. © = = = | 3 3 4 E a 
2-/2|)2 | 2 |25|.5/]2-] 2, | Ze | Ze] 8 
N a 5 ad a} a Oo & S wn 
Ba | g z € |aAe/S5]/8a] 425 | Sa | aa) & 
aes = a ac rallied 5 ‘So a oS =o & 
oat — . N — 
5 = 8 ~t 5 a | 8 3 3 u 
ay y S & x + 1S) oS ey S 
a) roa) ty 
Table 3 
Costs for the worker's salary for performing the operation 
3on = 1,872 Rub. (32 %) 
Tum.x. (29 %) T;, (3 %) 
Tum, (25%) Tn. (4 %) = 
Ton. (22 %) Tose. (2 %) Tomo(1%) Tn3.1 (1.5 %) T 3.2 ad %) Tip. (0.5 %) — 
Tex. o6c T, ope.o6c 
See Table 9 — See Table 8 See Table 9 — — 
(1.0 %) (1.0 %) 
Table 4 
Costs for manufacturing the product 
Ton. (22 %) 
To. Tinex.obe- T, 2.00+ 
Ts. (8 %) aie 
(14 %) (1 %) (2 %) 
Tusm- 
_ Tyem. (4 % Tynp. 3% i = — 
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— eee ae nae Ss oe og FM |e) a ze = 
HRYVl| NY] eS w S} 4s S BS) IN S| geo Se 
Table 5 
Costs for organizational preparation 
Tns.1 1.5 %) 
-©&] .€ | 8] .8] 28) s&] e&] g&] oF] g& | tH ge] es 
Za Za 4 =f SH fo 24 Rl eel 24 $4 54 fou 
“sl"el/Fs/"s/es]SsefFslfelss]ssgfeqeel] ig 
= i 
= 
5 Table 6 
a) : : 
av Costs for setting up the machine 
| 
$ Tn3.2 1.0 %) 
> os 8 . iio —~ = M1 rox — x om 
& s&] eel §8] FS sl] ese] ge ~& ye Se] gS] Re] Ss 
= sa | bal ent] eal So] Et] eS esl] eS) ea] $4) ES] §a 
“e|]ee ge dee eeles ("s/s liel/eg|/s| es 


162 


Bochkaryov PYu, et al. Comprehensive Assessment of the Manufacturability of Products 


ANNA Z 


WAY 

iN ¥. 
= AND 
al 


aw 


| 
1 
J 
t 


SI 


a ae ea 


. a oe 3 
: 90 


N 


Nk 


We Bee eT 

eh 2 | | 

TE. 

000 Sofod Tosrporren” 
2 Copano 


isan 
AL 


fete] = Kopnyc 


rere | mam 20/1 FOCT 977-88 


HDi 


vind 


WN 


tj 
Ne 


* Paszmep Gas cnpabok. 
4,116, +T4/2 


3 ** Pasmepe odecney. ump. 


4 Moxpeimue L{9.xp. — a i g22| 11 

5 Ycnemam xa npoywocm 2udpabnuveckum Gobnewuem 15 Mg xa a oe, , 

ZepMeMUYHOcM nxebmamuyeckum GabneHuem 10Mia YmeyKu u nagexue Onn Cr 2 

dabnenus He GonycKaomca 4 = ‘canponMaus” 
2.Lapal 


Fig. 8. Examples of design drawings of machined parts (““GAZPROMMASH” LLC) (part 2) 
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The results of 


calculations 


of 


the 
the manufacturability of parts on individual structural elements of production costs (site No. 1, “GAZPROMMASH” 


degree 


of 


relative 


influence 


of quantitative 


indicators of 


LLC) are presented in Tables 7—10. 
Table 7 
Impact of manufacturability quantitative indicators To 
(basic (process) time for the manufacture or processing of a unit of product) 
Estimated alalalst+{nuleflrfilalelalale|]s Weight 
indicators al hen ee “4 a a a a a 4 a a nS indicators 
1.1 =< Bo]? 1 0 0 0 1 2 1 2 2 13 0.098485 
1.2 Oo;-| 1 0 0 2 0 0 1 1 2 2 9 0.068182 
1.3 Oo; 1] - 0 1 1 0 0 2 2 0 0 7 0.05303 
1.4 1); 24 2 — 1 1 1 1 1 1 1 1 13 0.098485 
1.5 2/2) 2 1 — 1 1 1 1 1 2 2 16 0.121212 
1.6 2/;,2) 1 1 0 = 0 0 1 1 0 0 8 0.060606 
1.7 2/0); 1 1 1 1 = 1 1 1 0 0 9 0.068182 
1.9 0; 2) 2 1 1 2 2 = 1 1 0 0 12 0.090909 
1.10 0; 2] 0 1 1 1 1 1 = 1 1 1 10 0.075758 
1.11 1/1) 0 1 1 1 1 1 1 — 1 1 10 0.075758 
1.12 0; 0} 2 1 2 2 2 1 1 1 1 13 0.098485 
1.13 0; 0} 2 1 0 2 2 2 1 1 1 — 12 0.090909 
Table 8 
Impact of quantitative indicators of manufacturability Tynp 
(time to supply tool to workpiece) 
Estimated i a on} + nl ol rn]jalelrz Q o| 3 Weight 
indicators a i ole nen hens nen eee ee indicators 
1.1 — 2 2 O;1}] 1} 2 )2 7] 2 2 2 2 18 0.136364 
1.2 0 = 0 0 ;0; 1 1/1 1 1 1 1 7 0.05303 
1.3 0 2 = 2/14 1 1/1 1 2 2 2 15 0.113636 
1.4 2 2 0 -— | 07] 1 1/1 1 1 0 0 9 0.068182 
1.5 1 2 1 2);-}|2);2 7] 1) 2 2 2 2 19 0.143939 
1.6 1 1 1 1;0};-/; 1) 1 1 1 1 1 10 0.075758 
1.7 0 1 1 1/;0} 1/;]-) 1 1 1 1 1 9 0.068182 
1.9 0 1 1 1 1} 1 Le) =). 2 2 2 2 14 0.106061 
1.10 0 1 1 1/;0); 1 1/0] - 1 2 2 10 0.075758 
1.11 0 1 0 1/;0] 1 1/0] 1 = 2 2 9 0.068182 
1.12 0 1 0 2/0) 1 1; 0] 0 0 = 1 6 0.045455 
1.13 0 1 0 2/0) 1 1/0] 0 0 1 - 6 0.045455 
Table 9 
Impact of quantitative indicators of manufacturability 7n.3./ 
(time limit for organizational preparation) 
= Weight 
Estimated indicators | =| TV] @] SV] 2] eS) 5] 2] a = = 4 a 5 indicators 
1.1 —- | 2 | 2 1 1 17}; 2];2);2);2)2)2)/ 2 21 0.135484 
1.2 Oo;-/;0)]1)0/;/;0;,0;,0;0;,0;)0)7 07 0 1 0.006452 
1.3 Oo; 2)-) 1 1 1 1 1 1 1 1]; 0] 0 10 0.064516 
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= Weight 
Estimated indicators | = | @ | © = | De gd es ie le = zs 4 a 5 indicators 
1.4 1 1 1 0;/0/0);0;0;0;]0;] 07) 0 3 0.019355 
1:5 1 |} 2 1 O;,H 1 2),2),2) 2/2 )2 12 ) 2 20 0.129032 
1.6 1 |} 2 1 2/;0/;-];,2;);,2)2)2)/2) 1 1 18 0.116129 
1.7 0 | 2 1 2/0/]0j,- 4) 1 1 1 1 1 1 11 0.070968 
1.8 0 | 2 1 2/0/0)1;)-//)] 1 1 1 1 1 11 0.070968 
1.9 0 | 2 1 2/0/07] 1 -— | 1 1 1 1 11 0.070968 
1.10 0 | 2 1 2/0/07] 1 lL} = |} 1 1 1 11 0.070968 
1.11 0 | 2 1 2/0/07 1 1 1}-/ 1 1 11 0.070968 
1.12 0; 2/2 ;/2)04, 1 1 1 1 1; -) 1 13 0.083871 
1.13 O;2/;/2;/2 7) 1 1 1 1 1 1 1] - 14 0.090323 
Table 10 
Impact of quantitative indicators of manufacturability Acm 
(costs for the use of process facilities) 
- Weight 

Estimated indicators *% a = 5 4 a 5 indicators 

1.8 = 1 1 2 1 1 6 0.2 

1.9 1 = 2 1 1 1 6 0.2 
1.10 1 0 0 0 0 1 0.033333 
1.11 0 1 2 = 1 1 5 0.166667 

1.12 1 1 2 1 — 1 6 0.2 

1.13 1 1 2 1 1 = 6 0.2 


1.1 — material machinability index; 


1.2 — part design complexity index; 


1.3 — coefficient of accuracy and surface roughness of the part; 


1.4 — indicator of unification of structural elements; 


1.5 — material usage rate; 


1.6 — indicator of the possibility of manufacturing a given range of parts; 


1.7 — indicator of the use of production system capabilities; 


1.8 — indicator of the manufacturability of the part by the uniformity of process facilities; 


1.9 — indicator of predicting the level of loading of process facilities when processing a given range of parts; 


1.10 — indicator of multivariate decision-making when designing a process; 

1.11 — indicator of multivariate decision-making in the implementation of processes; 

1.12 — indicator of the manufacturability of the part, reflecting the possibility of observing the principle of unity of 
bases under the process development in terms of the surface of the part that is the main design base; 


1.13 — indicator of the manufacturability of the part, reflecting the possibility of observing the principle of unity of 


bases in the development of the technological process in terms of the surfaces of parts that are auxiliary design bases. 


Based on the presented models and previously known dependences of the calculation of quantitative indicators, an 
assessment of the manufacturability of parts was carried out. At the same time, the software developed and registered by 
the authors was used. Thus, taking into account the information about the real state of the production system, the 


adjustment of design documentation, range, sequence of implementation of the manufacture of individual groups of parts 


and production planning was performed. A comparative analysis of the calculation results is presented in Table 11. 
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Table 11 
Comparison results 
No. 1 option 2 option Performance 
Total time for the manufacture of products of 20 
1 items with an annual production program of 61193.42h 53073.35 h 15 % 
16,600 pcs. 
> Estimated number of equipment participating in 23 V7 35 % 
the process 
3 Operation factor 0.67 0.72 71% 


Discussion and Conclusion. The results of the presented theoretical studies and their approbation under real 
production conditions allowed us to propose a method for assessing the manufacturability of parts. It provides for a 
comprehensive assessment based on the developed analytical dependences for determining the weighting 
coefficients that characterize the significance of each indicator of manufacturability from the standpoint of the 
efficiency of the machining system. A distinctive feature and scientific novelty of the work is the consideration of 
the actual emerging production situation when assessing manufacturability. This makes it possible to use this design 
procedure not only traditionally at the initial stages of production planning, but also at the stages of implementation 
of processes for the purpose of rational organization of the production. 

The developed formalized models create the basis for complete sequential automation of design actions when 
evaluating the manufacturability of products, and provide prerequisites for constructing a promising intelligent 
system of predicting the efficiency of manufacturing parts in a particular production and making informed 


organizational and engineering decisions. 
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3aneneHHoll 6K1a0: 

IL.}O. boukapep — Hay4Hoe pyKOBOCTBO, pa3paboTKa MeTOAH4ecKOrO TOsxXosa OWeHKH BUMAHHA 
KOJIM4CCTBCHHBIX NOKa3aTesel TEXHOJOFMYHOCTH, aHalu3 pe3yIbTAaTOB UCcIeqOBaHHU, WOpaboTkKa TeKcTa. 

P.J]. Kopones — pa3paboTKa Moyelei B3aMMOCBaA3H WoKa3aTesel TeXHOMOrHYHOCTH HU 3aTpaT 
IIPOH3BOJ,CTBCHHOTO BPeMeHH, MpOMBIMIeHHad ampobamua H OOpaboTKa pe3yIbTATOB IKCHeCPHMeCHTOB, MOArO TOBKa 
TeKcTa. 


JLT. Boxopa — onpeyenenue coctaBa HM OleHKa NOKa3zaTelei NpOM3BOACTBeHHOM TeXHOJIOrH4HOCTH. 
Kond®aukm unmepecoe: aBtTopbl 3aABIAIOT OO OTCYTCTBHH KOHDJINKTa HHTepecos. 


Bce aemopbl npowmalu u odo6puanu OKOHYaMebHvIU 6apuaHmM PyKONUCU. 
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Abstract 

Introduction. The prediction module generates possible future trajectories of dynamic objects that enables a self-driving 
vehicle to move safely on public roads. However, all modern prediction methods evaluate their performance only under 
urban conditions and do not consider their applicability to the domain of rural roads. This work examined the adaptability 
of existing methods to work under rural unstructured conditions and suggested a new, improved approach. 

Materials and Methods. As a solution, we propose to use a neural network that includes the following submodules: a 
graph-based scene encoder, a multimodal trajectory decoder, and a trajectory filtering module. Another proposed feature 
is to use an adapted loss function that penalizes the network for generating trajectories that go beyond the drivable area. 
These elements use standard practices for solving the prediction problem and adapting it to the domain of rural roads. 
Results. The presented analysis described the basic features of the prediction module in the rural road domain, showed a 
comparison of popular models, and discussed its applicability to new conditions. The paper describes the new approach 
that is more adaptive to the considered domain of study. A simulation of the new domain was performed by modifying 
existing public datasets. 

Discussion and Conclusion. Comparison to other popular methods has shown that the proposed approach provides more 
accurate results. The disadvantages of the proposed approach were also identified and possible solutions were described. 


Keywords: trajectory prediction, behavior prediction, neural networks, self-driving cars, artificial intelligence, 


autonomous cars 
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AHHOTalna 

Beedenue. buaroyapa MoxyI10 MpornosMpoBaHuA TpaeKTOpHii ABWKeHHA THHAaMHYeCKUX OOBeKTOB OeciHJOTHBIM 
aBTOMOOMJIb CllocobeH Oe3s0MacHo ABHraTbedA 10 WoporaM oOulero Mob30BaHHA. OWHAKO BCe COBPeMeHHbIe MeTObI 
TIpOrHO3HpOBaHHA OLWCHHBAaIOT MPOW3BOUTeIbHOCTh TOJIbKO B TOPOCKHX YCJOBHAX MH He PaCCMaTPHBalOT CBOIO 
IIPHMe€HUMOCTb K JOMeHy Mpoceso4HBIx AZopor. Len DZaHHoro UcceqOBaHHA 3aKkOUAeTCA B AHaIM3e alalTHBHOCTH 
CYUIeCTBYIOWJMX MeTONOB MIpOrHo3HpoBaHHA MU pa3paOoTKe MOAXoya, KOTOpbIM OyeT AeEMOHCTpHpoBaTh JIy4LLy10 
IIPOH3BOHTEJILHOCT IIpH paOoTe B HOBBIX YCJIOBHAX. 

Mamepuanoi u memoooi. B kauectBe pelieHua MpeslaraeTcA MCMONb30BaTb HEHPOHHYIO CeTb, BKJIOYAIOLIyIO B cebs 
crelyroulHe MOAMOAyIM: rpapoBblit KOAMPOBLIHK CI[CHbI, MYJIbTHMOMaIbHbIM JeKOAMpOBWIHK TpaeKTOPpHi, MOJyJIb 
(puIbTpalyun Tpaektopui. TawKe mpesiaraeTca NPHMeHUTb alanTHupoBaHHylo PYHKUMIO MOTepb, KOTOpad WITpadpyer 
CeTb 3a TeHepallMiO TpacKTOpH, BEIXOAAMUX 3a TpaHHUbl AOpoxkHoro MowoTHa. JJaHHble IIeMeHTbI 3aeHCTByIOT 
paciipocTpaHéHHble MpaKTHKH PeWIeHHA 3aa4n NPOTrHO3MPOBAHHA, a TakOKe aflalITHpyIOT eé WI WOMeHa MIpOcesO YHBIX 
Wopor. 

Pezynsmamot uccredosanua. ([poananu3upoBaHbl OCHOBHBIe OTIIMYHA H YCIOBHA paOOThI MOJYIA IpOrHOSHpOBaHHA B 
yCJIOBHAX MIPpOCeIOUHBIX Jopor. BenosHeHa CHMYJIAUHA HOBOTO JOMeHa IlyTEM MOAMUuKAalHH CyLecTBYIOWIMX HabopoB 
aHHbix. IpopexeHo cpaBHeHHe MOMMyIAPHbIX MeTOJOB MporHo3sHpoBaHHA HM OWeHeHa HX IPHMCHHMOCTb K HOBbIM 
ycnoBusaM. IIpezcTaBiieH HOBbIM, OosIee alanTHBHbIM K HOBOMY JOMeHYy, MOAXod. 

O6cyscdenue u 3akimio4uenue. IIpopeqeHHoe cpaBHeHve C JIpyrHMM MOMyJIApHbIMH MeTOJAaMM TOKa3bIBaeT, 4TO 
Tipesox*KeHHOe aBTOpaMH pelieHHe obecreauBaeT Oosee TOUHbIC pe3yIbTaTbI MporHo3sMpoBaHHna. TawKe Obi 


BBIABJICHbI HELOCTATKH MpeAIO#xKCHHOTO MOAXOa HW OMMCaHbI BO3MO2KHBIC IlyTH HX YCTpaHeHHsA. 


KoroueBbie cJI0Ba: TIpOrHOSHpOBaHHe TpaeKTopui, TIPOrHO3HpOBaHHe NOBeTCHHA, HelpoHHble ceTH, OeciHHIOTHBIC 


aBTOMOOMJIN, MCKYCCTBeHHbIi UHTCJVICKT, ABTOHOMHBIC aBTOMOOMIM 


BuaarogapHocrn: aBTOPbI BbIpaxKaroT OlaroyapHOCTb WeHTpy OecHHJIOTHBIX TEXHOJIOMM yHHBepcuteta VWnHonomne 3a 


TIOMOIMb B HpOBeTeCHHH HCCIeTOBaHHA. 


Aaa untTupopanna. Veaxos C.A., Pamuy b. IIporno3supopanue noBpeyeHua y4acTHHKOB JOPOXKHOTO BYWKeHHA B 
YCIOBHAX IpOCeNOYHBIX WOpor Ayia OeciMMOTHBIX aBTOMOOUIeH. Advanced Engineering Research (Rostov-on-Don). 
2023;23(2):169-179. https://doi.org/10.23947/2687-1653-2023-23-2-169-179 


Introduction. The latest achievements in the field of artificial intelligence (AJ) are being actively implemented in 
various areas of activity. One of such achievements is autonomous vehicles (AV). The current research was aimed at 
creating algorithms that allow AV to move safely on public roads. This will significantly reduce the number of road 
accidents [1]. 

The scientific community has already identified the basic modules of an autonomous vehicle. One of them is a 
system for predicting the future behavior of road users (agents) [2]. A clear understanding of how the environment will 
develop and in which direction dynamic objects (pedestrians, cars, cyclists) will move is urgently needed for AV to 
search and use a safe and effective trajectory of movement. 

Numerous scientific papers are devoted to the problem of predicting such trajectories [3-12]. However, there is 
currently no active research on the application of existing methods outside of urban conditions. And this is extremely 


important, since autonomous cars will be used on country roads, too [13]. Urban conditions are highly structured: cars 


Ivanov SA, et al. Predicting the Behavior of Road Users in Rural Areas for Self-Driving Cars 


mostly follow traffic lanes, and pedestrians move through special zones. In this sense, the area of country roads is the 
complete opposite, which means that it will have additional difficulties during development. In the given paper, attention 
is focused on these difficulties: the existing predicting methods and their applicability to new, less structured conditions 
are considered. 

The objective of the study involved: 

—analysis of the major differences and working conditions of the prediction module under the conditions of country roads; 

— simulation of a less structured country road domain by modifying the existing datasets; 

— comparison of modern prediction methods, including their applicability to new conditions; 

— description of the new approach and proof of its higher accuracy in comparison to other prediction methods. 

Materials and Methods. At first glance, it would seem that the domain presented is a simpler version of urban 
conditions due to the fact that country roads are characterized by less traffic. However, the absence of complex multi- 
level junctions, special traffic-free zones, a large number of signs, markings, etc., makes the domain of country roads less 
structured, i.e., fewer rules and specific traffic patterns increases the randomness and reduces the predictability of the 
behavior of cars and pedestrians. 

The following features of the country road domain influence strongly on the selection of the architecture of the 
prediction module: 

— crossroads. Undoubtedly, they are more simple on country roads in comparison to urban ones, but at the same time 
this fact of simplicity means that the model must take into account multimodality and assess the probability of choosing 
each possible direction of movement at the crossroad when the agent approaches it; 

— country roads do not have lane markings, pedestrian crossings, bike paths, etc. Instead, the HD map will contain 
only information about the boundaries of the roadway. Therefore, the stage of encoding the scene should take into account 
this feature to describe the surrounding context more effectively; 

— pedestrians and cyclists will move along the same road with conventional and autonomous vehicles. Therefore, the 
model should be adaptive for predicting the future trajectories of both cars and pedestrians/cyclists. 

The prediction module implies the presence of AV recognition, tracking and localization systems and their accurate 
operation. The authors of the article use the Argoverse dataset, which stores the required records of the operation of all 
systems in a convenient form [14]. 

The dataset consists of recordings of road scenes observed on the streets of Miami and Pittsburgh, USA. Each of the 
entries contains a local part of the terrain map (lane boundaries, roads, pedestrian crossings) and a list of all recognized 
agents, including the current position and movement history of each of them. Each of the records is divided into two parts: 
two seconds of the observation history and three subsequent seconds for which prediction is made (prediction horizon). 
Data on the future movement of objects is also available and used to calculate the accuracy of prediction methods and 
model training. 

Information about agents is presented in a discrete format. The time interval between measurements is fixed, in this 
work it is equal to 0.1 seconds (10 Hz). 

For each moment of time ¢, the module receives the observation history or for each detected agent i. The observation 
history consists of the agent's current and past states, where each of the states s, is a 2D position in the global coordinate 
system. The authors make the assumption that the height information is redundant. 

The dataset also provides access to an HD map that contains information about road borders and roadway, pedestrian 
crossings. To simulate the domain of country roads, the dataset was modified in such a way as to exclude all information 
from road maps, except for the boundaries of the roadway D. This reduces the amount of information about the road 
context and complicates the task of prediction. 


Hence, the context of the scene is represented as 
CHATS? BP SPD), (1) 


where k — the total number of tracked agents on the scene. 
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This approach implies predicting the trajectory for only one agent per execution, therefore further S,’ is treated as 


S” for simplification. To generalize the model for all recognized agents, it is required to repeat the proposed approach 
for all k agents on the scene. The agent for which the prediction is currently being made is considered a target agent. 


To assess the accuracy of prediction methods, the dataset contains recorded future trajectories S/ for each target agent: 
SS i Sina aiel (2) 


where H indicates the number of next time steps. In this case, parameter H will be equal to 30, since the planning horizon 


is three seconds with a sampling frequency of 10 Hz. 

The domain of the prediction module is multimodal, i.e., the future behavior of agents may differ significantly in 
absolutely identical traffic situations. Let us say, a car approaching a crossroad may continue straight ahead or make a 
turn. To take this into account, it is required to generate M possible future trajectories and M probabilities of each of them 
at the model output. 


Therefore, the purpose of the prediction module is to create function f, that takes the context of the scene c as input 


and generates M pairs of possible future trajectories and their probabilities: 
FC) = SF Si o9St PaDae sat (3) 


Here, at least one generated trajectory S/ should be as close as possible to the real trajectory S/ , and the probability 
of its execution p should be close to unity. 
Model architecture. The proposed approach involves the use of a neural network consisting of submodules of scene 


encoding, decoding and filtering trajectories. The architecture of the system is shown in Figure 1. 
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Fig. 1. System architecture 


A neural network adapted to the new conditions, based on a vector representation, is responsible for encoding the 
scene. This selection is due to the fact that on country roads, the HD map will contain a limited amount of information 
(only the roadway boundaries and the history of observations of dynamic objects). Popular methods represent the context 
of a road scene c in an image format and process it using convolutional neural networks. However, vector coding avoids 
the overhead associated with image generation [4—5]. 

The presented encoder is based on the VectorNet model, but its input data format has been modified to receive 
information only about the boundaries of the roadway and the state of the agents. [3]. This encoder represents the 
boundaries of the road and the state of agents using polylines, which are further processed by a graph neural network. 
This provides encoding the interaction between polylines. Details of the implementation are described in paper [3]. 

A trajectory decoder is a task of regressing several possible trajectories and generating a set of probabilities. To solve 
this problem, a multilayer perceptron model is used. The decoder implementation is inspired by the MTP model [4], 
however, the authors of the article propose a different formula for calculating the best trajectory m* from the set of 
M trajectories. It is also proposed to use an additional mechanism that penalizes the model for predictions that go beyond 


the area of movement. 
The authors of the original MTP model propose to train a multilayer perceptron using the loss function that represents 


the sum £,,, and £,,,,,, where: 


lass ? 


b= ES" 85) () 


Teg 
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M 
Ltass = -yi log Pn (5) 
m=1 


In this case, £,,, —- mean-square error between the real trajectory S ! and the best trajectory m* of M generated. 


e(s'.s/)=2>(s-s,) » ©) 


i=l 
where s;— the agent's actual future position at time i, and 5 — the predicted future state of the best trajectory m’*. 


L, 


asx ~~ 1088 function based on cross entropy, which increases the probability of executing the best of the predicted 
trajectories m* to 1 and reduces the probability of other trajectories to 0. 

[, is a binary indicator equal to 1, if condition c is true, and 0 — otherwise. 

In the original article, the best of the predicted trajectories m* is defined as the one that has the minimum value of the 


root-mean-square error in comparison to the real trajectory: 


m=argmin (5? 82). (7) 


m 
m 


The authors of the article suggest using the following modification: 


m 


m=argmin (sts!) (8) 


meA 


where A — subset of generated trajectories that has a similar final direction to the real trajectory S/ . 

The idea is to remove from consideration trajectories in which the final direction of the agent differs significantly 
from the direction in the real trajectory when calculating the best trajectory m*. If the difference in directions features 
less than certain threshold y, then the generated trajectory is considered correct, i.e., mA. In the case under consideration 
y=30°. Therefore, the best trajectory m* should have a similar final direction and the lowest value of the loss function. 

This work also involves prior knowledge of the domain to achieve greater convergence of the model [15]. Since only 
information about the roadway boundaries is available from the HD map when driving in the domain of country roads, 


an additional variable is introduced — £,, into the loss function. Thanks to it, the model will penalize the predicted 


trajectories that go beyond the road in cases where at least one state is s; € D. The model penalizes only the best trajectory, 


since only in this case, it is possible to determine the direction of error reduction by approximating the best of the generated 
trajectories m* to the real trajectory S’. 


Thus, £,, is defined as: 


12 + \? 
Lin == YI.-(5-51) ’ (9) 
where /, is equal to 1, if s, ¢ D , and 0 — otherwise. 


The final loss function is defined as 


L= Lig tO Lo +B Lips (10) 


lass reg 
where a and B — neural network hyperparameters used for training. In this case, both of these parameters are equal to 0.5. 

To filter similar and duplicate trajectories, the proposed approach uses the filtering of a finite set of trajectories M at 
the final stage. This module is required because in some cases, the number of possible agent trajectories may be less 
than M, e.g., when a car is moving along a straight road at a constant speed, the model can generate only one trajectory: 
the car continues to move straight. However, the need to generate exactly M trajectories will result in the situation when 
all predictions are similar to each other. 

The proposed filtering is based on the final direction and positions of states s;: if the direction and the sum of the 
deviations between states s; of the real and generated trajectories are less than the threshold value o, then the trajectories 
are considered similar. The authors average each state of the trajectories and sum up the probabilities of the trajectories p;. 

This approach was implemented in the Python programming language on the PyTorch deep learning framework. The 
model was trained on GeForce RTX 2080 Ti graphics card for 40 epochs, the training took four hours. 
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Research Results. To assess the accuracy of prediction models, this section applies widely used metrics for the 
trajectory prediction problem: average displacement error, ADE, final displacement error, FDE [6], MissRate (MR), and 
Offroad rate (OR). 

For multimodal cases with the generation of several trajectories, ADE and FDE are taken as the minimum ADE and 
FDE among M trajectories (the trajectory with the lowest metric value) [5]. 

The prediction is considered “missed” if ADE metric of the generated trajectory is more than two meters. OR metric 


is calculated as the percentage of trajectories in which at least one state s; goes beyond the range of motion D. 
To visualize the context of scene c, as well as the real and predicted future trajectories Sand S’ a script in the Python 


programming language was implemented using the Matplotlib library. 

This section compared the operation of several different methods in the case of an unstructured domain. The following 
methods were used in comparison: 

— Kalman filter; 

— Single trajectory output — the proposed scene encoder with the generation of a single trajectory; 

— Fixed set classification — the proposed scene encoder with the reduction of the task to classification among 
predefined trajectories: by sets of 64 and 415 predefined trajectories; 

— Proposed approach. 

Table 1 presents comparison of the accuracy of the methods when working under unstructured conditions. Several 


methods are compared, including the proposed approach. 


Table 1 
Comparison of models in the unstructured domain of work 

Method Modes | ADE, FDE, ADEs FDE, MR2, MR2. OR 
Kalman filter 1 3.78 8.05 3.78 8.05 0.89 0.89 5.89 
Single trajectory output 1 3.12 6.75 3.12 6.75 0.89 0.89 3.26 
Fixed set classification 415 3.27 7.00 1.74 3:57 0.84 0.52 3.61 
Fixed set classification 64 2.6 5.63 1.52 2.91 0.82 0.49 2.58 
Proposed approach 6 2.36 5.29 1.32 2.55 0.78 0.38 1.84 


Kalman filter. The simplest way to predict behavior is to obtain the current state of the object (current lane, speed, 
direction, etc.) and extend this state to future steps based on some assumptions, e.g., that the car will continue to follow 
its lane or will have a constant speed and/or acceleration. Another popular method for such tasks is to use the Kalman filter [12]. 

According to Table 1, the Kalman filter works worse than all the presented methods based on neural networks. 

Figure 2 shows two cases. In the first case, the Kalman filter successfully performs prediction because the vehicle is 
moving straight, without any turns or speed variation. In the second case, the Kalman filter mispredicts due to lack of 


knowledge about the context of the traffic situation. 
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Fig. 2. Example of predictions using the Kalman filter. Dotted lines — roadway boundaries, 


red lines — target agent with history of observations, blue — other agents, green — real trajectory, yellow — predicted 


trajectory, red crosses indicate predicted states outside the roadway 
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Single trajectory output. This method involves the use of a graph scene encoder, which is identical to the one used in 
the proposed approach. The output of the network implies the generation of only one trajectory. This model is trained 
using the root-mean-square loss function. 

As shown in Table 1, the neural network, even with the generation of a single trajectory, demonstrates better results 
in comparison to the Kalman filter. 

Figure 3 shows the visualization of this prediction method operation. The image on the left shows that the model can 
successfully predict the agent's turn. The image on the right shows that generating one trajectory is not enough. The neural 
network tries to imagine both possible outcomes: going straight and turning right. As a result, the model outputs the 


average of the two outcomes. 
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Fig. 3. Example of generating a single trajectory. Red line shows the target agent with 


history of observations, green — real trajectory of movement, yellow — predicted trajectory 


Fixed set classification. The implementation was inspired by the CoverNet prediction method [5]. This model consists 
of a proposed vector scene encoder, followed by a different trajectory decoder. The decoder is a classification task based 
on a predefined set of trajectories consisting of physically realizable vehicle trajectories with sufficient coverage. Two 
sets were created for experiments: of 415 and 64 possible trajectories. The second set has the same coverage as the first, 
but provides a lower density of trajectories. Detailed information about the sets of trajectories is contained in paper [5]. 

The visualization of the work is shown in Figure 4. The classification model successfully copes with multimodality at 


crossroads, but in some cases, the lack of sufficient coverage by a set of trajectories negatively affects the results. 
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Fig. 4. Example of prediction using a classification model. Red lines represent the target agent < 

with history of observations, green lines — real trajectory. M predicted trajectories a 

with different probability of execution pi are presented using red-yellow hues x 

| 

As shown in Table 1, this method works more accurately than generating a single trajectory, but worse than the Be 
proposed approach. In addition, increasing the density of the set of trajectories by using a set of 415 trajectories did not E 
improve the results. The authors attribute this to the presence of noise in the dataset, which comes from the tracking = 


system used in the data collection. 
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Proposed approach. The proposed approach eliminates the disadvantages of all the methods described above. This is 
a multimodal forecasting method that does not suffer from the limitations of a predefined set of trajectories. 
Moreover, according to Table 1, the proposed approach surpasses all other methods in all indicators. As shown in 


Figure 5, the method successfully captures two possible outcomes at the crossroad: driving straight or making a turn. 
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Fig. 5. Example of prediction using a classification model. Red lines represent the target agent 
with history of observations, green lines — real trajectory. M predicted trajectories 


with different probability of execution p; are presented using red-yellow hues 


Figure 6 shows an example of filtering similar trajectories in the case of a single possible outcome. The probability 
that the agent will complete the initiated turn is close to 1, since he is already in the process of turning. Therefore, in this 


case, the probability of other outcomes is close to 0. The proposed module successfully filters similar trajectories. 
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Fig. 6. Filtering effect. The entire set of predictions is shown on the left, and only filtered set — on the right 


Limitations. Although the authors of the original article on the MTP model [4] indicate that their method solves 
the problem of mode collapse, the experiments conducted by the authors of this article do not confirm this. The 
problem still occurs in some cases. It is assumed to be due to the following features: the loss function does not 
penalize the neural network for generating all possible trajectories that the target agent can execute, as long as the 
best of them is as close as possible to the real trajectory. But also, the model does not encourage the network in any 
way to predict a variety of possible trajectories. Therefore, it is advantageous for the network to make several similar 
predictions in one direction, in which it is more confident than to make one prediction for each possible trajectory. 

One of the possible solutions to this problem may be the use of a trajectory decoder presented in the TnT, 
DenseTnT models [10—11], which imply the generation of final goals at the first stages of work. In these models, 
all possible final goals for the agent are generated first, and then trajectories that describe the movement from the 
starting position to each of the goals, are generated. This provides filtering out similar final goals in the early stages, 


and preventing the mode collapse. 
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Discussion and Conclusion. In the work performed, modern methods of solving the trajectory prediction 
problem are investigated. The adaptability of the methods to unstructured road conditions — country roads, is 
considered. Insufficient accuracy of the methods is established, and a new approach to predicting is proposed. 

The proposed approach is based on the VectorNet and MTP models, but has been adapted for the country road 
domain. In addition, a trajectory filtering module and an additional mechanism for the loss function, which penalizes 
trajectories for going beyond the movement zone, are proposed. 

The presented comparison shows that the proposed approach is superior to other popular methods. 

Limitations of the MTP approach have been identified: the output data still tends to mode collapse. The 
suggestion for further modifications is to use methods that generate the final goal at the early stages of prediction 


and thus are less susceptible to regime collapse. 
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Abstract 

Introduction. The 2020s were marked by the emergence of a new generation of computer simulators using augmented 
reality. One of the promising advantages of augmented reality technology is the ability to safely simulate hazardous 
situations real-world. A prerequisite for realizing this advantage is to provide the visual coherence of augmented reality 
scenes: virtual objects must be indistinguishable from real ones. All IT leaders consider augmented reality as a next 
“big wave”; thus, the visual coherence is becoming a key issue for IT in general. However, it is in aerospace 
applications that the visual coherence has already acquired practical significance. An example is Boeing's development 
of an augmented reality flight simulator, which began in 2022. Visual coherence is a complex problem, one of the 
aspects of which is to provide the correct overall coloration of virtual objects in an augmented reality scene. The 
objective of the research was to develop a new method of such tinting. 

Materials and Methods. The developed method (called spectral transplantation) uses two-dimensional spectral image 
transformations. 

Results. A spectral transplantation technology is proposed that provides direct transfer of color, brightness, and contrast 
characteristics from the real background to virtual objects. An algorithm for automatic selection of the optimal type of 
spectral transformation has been developed. 

Discussion and Conclusion. Being a fully automatic process without recording lighting conditions, spectral 
transplantation solves a number of complex problems of visual coherence. Spectral transplantation can be a valuable 


addition to other methods of providing visual coherence. 
Keywords: computer simulators, augmented reality, visual coherence 
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AnHOTalnA 

Beedenue. 2020-e roqbl O3HAMCHOBAIMCh TOABJIGHHCM HOBOTO MOKOJICHHA KOMIIbIOTEPHbIX TPeHaKePOB C 
IIPHMeHeEHHeM TeXHOJIOrHM AOMOUHEHHOM peasIbHOCTH. OHO 43 NpPeHMYyecTB JaHHOM TeEXHOJOrHH — BO3MO2KHOCTb 
Oe30MacHOro MOJCIMpOBaHHA ONaCHbIX CHTyalMii B peasIbHOM Mupe. HeoOxogHMbIM yYCIOBHeM HCHOIb30BaHHA 3TOTO 
IIpeHMyllecTBa ABIAeTCA OOECHeYeHHe BU3YaIbHOM KOTePeCHTHOCTH CI[eH AOMOMHEHHOM peasIbHOCTH: BUPTyaJIbHble 
OObEKTEI OJDKHEI ObITh HEOTIIMYMMBI OT peasbHBEIx. Bce muposbie IT-nuyepbr paccMaTpHBaIoT JOMOJIHeHHYIO 
peaIbHOCTbh KaK CJI€MYIOWyIO BOJHY paMKasIbHbIX H3MeHeHH B WHppoBol cpefe, NOSTOMy BH3yasIbHaa 
KOTe€peHTHOCTb CTaHOBHTCA KJIIOUEBBIM BOIIpOcoM Wa Oyzyulero IT, a B aspPOKOCMHYECKHX MPHIOKCHUAX BH3YAaJIbHaAd 
KOrepeHTHOCTh yxe MpHoOpena upakTHyeckoe 3Ha4yeHHe. IIpuMepoM MOxeT CILyKHTb pa3spaOoTKa KopmopayHel 
Bouur MuuOTCKOrO TpeHaxepa C QONONHEHHOM peasbHoOcTbIO (2022). Bu3yanbHad KOrepeHTHOCTb — CJIO2%KHaA 
KOMIVIeKcHad mpoOsemMa, OJHHM H3 acieKTOB KOTOpOM sABUIAeTCA OOeCHeYeHHe KOPpeKTHOM kKoMOpucTuyecKoH 
TOHMPOBKH BUPTYaJIbHbIX OOBEKTOB B ClleHe AONOUHEHHON peambuoctTu. Lien padoTsr — pa3paboTKa HOBOrO MeTOa 
TaKOM TOHHPOBKH. 

Mamepuaavi u memoooi. B pa3spa0oTaHHOM MeToje (Ha3BaHHOM CII€KTpasIbHOM TpaHCiaHTalMen) UCMONb3YIOTCA 
J{BYMEPHBle CieKTpaJIbHble IpeoOpa30BaHuA H300paxKeHHH. 

Peszynomamoti uccredoeanua. I[peaioxwena TeEXHOOrMA CHeKTpabHOM TpaHciiaHTayun, oOechewnBalolsad IpAMyto 
Tlepeqauy xapakTepHCTHK WBeTa, APKOCTH HW KOHTpacta OT peabHoro (:boHa K BUPTYaJIbHbIM OOBeKTaM. Pa3spaboTaH 
ayIFOPUTM aBTOMATHYECKOTO BbIOOpa ONTHMAJIbHOTO BH Ja CIeKTpasIbHOrO MpeoOpa30BaHHA. 

O6cyacdenue u 3akniouenue. byyun MONHOCTbIO aBTOMaTH4eCKHM TpoweccoM 6e3 peructTpaluu ycroBui 
OCBeINeHHOCTH, CII€KTpasIbHad TpaHciaHTauuA pelllaeT pA CIOKHIX TpoOseM BH3yaIbHOM KOrepeHTHOCTH. 
CrekTpasIbHad TpaHCiaHTallHA MOKET CTaTb ICHHbIM JOMOJHeEHHeM K J{pyrMM MeTOJaM OOecrieyeHHA BH3yasIbHOH 


KOTepeHTHOCTH. 
Ksr0ueBble CJI0Ba: KOMIIBIOTCpHble TpCHaKepbl, WONOJHCHHAA PealsIBHOCTh, BH3yasIbHawA KOTCPpCHTHOCTh 


BuaarofapHoctu: aBTop BbIpaxaeT OnarogapHocTp A. Tepenuu (Inglobe Technologies Srl, Uexxano, Uranus) 3a 


TOWWepxKy B pa3paboTKe MmporpaMMHOro oOecrieyeHHaA. 


Aaa waTupopanna. CopOyxos A.J]. BusyanbHad KorepeHTHOCTh B JOMOMHeEHHOM peambHoctu. Advanced Engineering 
Research (Rostov-on-Don). 2023;23(2):180—190. https://doi.org/10.23947/2687-1653-2023-23-2-180-190 


Introduction. Modern simulators actually by default imply the use of virtual reality (VR). The advantages of this 
approach are well known; therefore, we will not dwell on them, but we will note a number of significant and, more 
importantly, insurmountable disadvantages due to the very nature of virtual reality technology. VR is a digital, discrete 
technology, while the real world is continuous. Therefore, modeling the real world in VR is inevitably associated with 
errors, which reduces the efficiency of training. However, for training systems, an even more serious negative aspect 
is that human decisions are largely based on subconscious consideration of numerous details of the real picture of the 
world. This process is fundamentally impossible to reproduce using purely computer technologies (e.g., VR) for two 
reasons: we still do not know (and are unlikely to ever know) what the mechanism of the human brain is. The latest 
speculations on the topic of artificial intelligence only confirm this. The details of the real world taken into account 
when making decisions are almost infinite in number, they arise randomly and are of quite a different nature (visual, 


acoustic, tactile ...). 
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The emergence of augmented reality (AR) training systems in the 2020s reduced the severity of this problematic 
situation. Examples are the development by Boeing of an augmented reality pilot simulator based on the well-known 
R6 ATARS project, which began in the fall of 2022, as well as a similar project launched by British BAE Systems or 
an air traffic control training simulator from this article. All the information wealth of the world around us in AR is 
presented explicitly and does not require modeling. But it is needed to solve the problem of visual coherence (VC) to 
realize the advantages of AR associated with the parallel presence of real and virtual objects in scenes: virtual objects 
must be indistinguishable from real ones. This article proposes a method for solving the problem of visual coherence 
in the framework of a project on the development of a training system for air traffic controllers. 

AR is a derivative form of VR. AR retains all the features of VR, but, in addition, as a hybrid technology, it has 
significant advantages arising from the parallel coexistence of virtual and real objects, which attracts the attention of 
developers to VC. Moreover, studies [1] show that among the negative psychophysiological consequences of using 
augmented reality devices, optical discomfort dominates, which occurs due to the difference in perception of real and 
virtual objects in the same scene due to the absence of VC. IT industry leaders see AR as the next “big wave” of 
revolutionary changes in digital electronics. Therefore, the VC problem is becoming a key one for IT as a whole, and 
these leaders show a growing interest in methods of solving it [2]. However, the problem of visual coherence has 
already acquired practical significance in aerospace applications. The authors encountered a VC problem when 
developing a training system for air traffic controllers: the rapid increase in the intensity of air traffic at airports caused 
an increase in the frequency of collisions of aircraft with other aircraft and airfield transport during ground 
maneuvering (>50 cases worldwide in 2018 before the outbreak of the pandemic). Air traffic controllers working on 
airport towers are not always ready to respond adequately to such emergency situations, which requires additional 
training. The most effective form of such training involves presenting the dispatcher with a situation of hazardous 
proximity of objects on the airfield, which is impossible with real objects, but can be absolutely safely implemented 
in augmented reality scenes. In our application, emergency situations were safely simulated using AR at a real airfield, 
while the virtual aircraft used should be indistinguishable from real ones. 

An exhaustive overview of the known VC methods can be found in [3]. According to the author's classification, all 
VC methods can be divided into two main classes: with the measurement of lighting parameters, and with the 
assessment of lighting conditions. In the first case, a mandatory procedure is a preliminary measurement of illumination 
conditions, carried out with the help of special equipment. This procedure is a long and labor-intensive process. It 
seems to be impossible if a pre-obtained image or video of the real world is used. In the second case, the complexity 
of reconstructing the lighting pattern from images causes assumptions and limitations, which makes the results 
ambiguous. Therefore, despite the impressive results obtained by researchers using the methods mentioned in the 
review [3], the VC level is still often insufficient, specifically, in AR scenes with real natural landscapes under ambient 
lighting conditions, which are typical for aviation applications. As the review of publications below shows, there is a 
shortage of research of this kind. 

This work was aimed at developing a universal and automatic method to provide direct transfer of color, brightness 
and contrast characteristics from a real background to virtual objects without digital 3D modeling, which was required 
in existing VC approaches. The method is based on the mathematical apparatus of two-dimensional spectral 
transformations, we called it “spectral transplantation”. 

The key results of this study are: 

— basic scheme for the spectral transplantation method, which provides a direct transfer of color, brightness and 
contrast characteristics from the real background to virtual objects. The method involves replacing a part of the 
spectrum of the image of the virtual world with the same part of the spectrum of the image of the real world, followed 
by an inverse transformation of the spectrum with the transplanted part; 

— algorithm for automatic selection of the optimal type of spectral transformation for use in spectral transplantation. 

It is important to note that VC depends on many factors: lighting, shadows, color tone, mutual reflections, surface 
texture, optical aberrations, convergence, accommodation, etc. Accordingly, various AR visualization techniques were 
used. In our case, VC is provided only for the factors of general illumination and coloring of virtual objects in AR. 
This is one of the VC challenges, especially for outdoor scenes. Therefore, spectral transplantation should be used in 


combination with other VC methods to achieve full VC. 


Gorbunov AL. Visual Coherence for Augmented Reality 


The list of sources in [3] includes 175 positions; this review includes almost all approaches to achievements in VC 
(with the exception of the latter, based on neural networks, discussed below). Therefore, here we will briefly describe 
some characteristic examples that correspond to the mentioned basic classes. 

Measurement of lighting conditions 

Using a light probe with diffuse bands between mirror spherical quadrants, P. Debevec and others [4] demonstrated 
how the full dynamic color range of a scene could be reconstructed from a single exposure. Based on the image obtained 
with the probe, the intensity of several light sources could be estimated by solving a simple linear system of equations. 
The results were used to render a virtual diffuse sphere. 

A. Alhakamy and M. Tuceryan [5] estimated the direction of incident light (direct illumination) of a real scene 
using computer vision techniques with a 360° camera attached to an AR device. The system simulated the light 
reflected from surfaces when rendering virtual objects. Then, the shadow parameters for each virtual object were 
determined. 

Assessment of lighting conditions 

S.B. Knorr and D. Kurtz [6] proposed a scheme for assessing lighting conditions in the real world based on a photo 
of a human face. The method was based on training a model of the type of face based on a database of faces with 
known lighting. The authors then reconstructed the most plausible lighting conditions in the real world in the basis of 
spherical harmonics for the captured face. 

We should mention work [7], which described a combination of measurement and evaluation of illumination. The 
authors measured the reflective properties of real objects using depth maps and color images of a rotating object on a 
turntable using an RGB-D camera. The shape of the object was reconstructed through integrating images of the depth 
of the object obtained from different viewpoints. The reflectivity of an object was determined by evaluating the 
parameters of the reflection model from reconstructed images of shape and color. 

The closest analogues of the proposed method are approaches that, like spectral transplantation, do not involve 
preliminary measurements of lighting and simulation of lighting conditions, scene geometry, surface reflection, and 
also provide for automatic processing. 

Among such analogues, there are methods of color transfer from image to image. Paper [8] presented a method for 
automatic transferring color statistics (averages and standard deviations) from the reference image to the tar get image. 
Additional parameters were used to avoid manual processing, which was required to determine the features of color 
transmission in cases where images had a strong difference in the color palette. These additional parameters combined 
the variances of the reference and target images. The authors of the article claimed that, although manual modification 
of these parameters was extremely rare, it was nevertheless sometimes necessary. In addition, the statistical nature of 
the method raised questions about the type and scope of statistics. Also, the ability of the method to process certain 
types of images (containing shiny objects, shadows) was not obvious. 

Xuezhong Xiao and Lizhuang Ma [9] presented an algorithm to solve the problem of color transmission reliability 
in terms of scene details and colors. The authors considered the preservation of the color gradient as a necessary 
condition for the authenticity of the scene. They formulated the problem of color transfer as an optimization problem 
and solved it in two stages — histogram matching and gradient-preserving optimization. A metric was proposed for 
an objective assessment of the efficiency of color transfer algorithms based on examples. 

The advantages of the developed method, in comparison to [8, 9] and their numerous analogues, are its versatility, 
fully automatic nature, and the ability to transfer not only color, but also all the main characteristics of the image using 
one simple procedure. 

The proposed method uses two-dimensional spectral transformations. Various types of images are optimally 
described by different types of spectral transformations (“‘optimally” — in the sense of matching visual perception for 
real and virtual objects). Actively used in digital image processing since the advent of digital television are the Discrete 
Fourier Transform, Discrete Cosine Transform, Hadamar Transform, S-Transform, and Karhunen-Loeve Transform. 

Materials and Methods. The scheme of the spectral transplantation method (the version using the Fourier 
transform [10]) is shown in Figure |. Frames of the real world (world frame — WF) and virtual world (virtual frame — 


VF) are used as input data (Fig. 2). This is natural for AR “video” (when the real world is observed through a video 
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camera). For “transparent” AR, when the real world is observed through transparent glasses, real images are captured 


using cameras located on AR glasses. 
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Fig. 1. Scheme of the spectral transplantation method. 


A version using the Fourier transform 


Fig. 2. Real world (WF) and virtual world (VF) video frames: 
a — WF, airport, cloudy weather; b — WF, airport, sunny weather; c — VF, 


virtual airplane. WF are small fragments (<25%) of images published on websites sydneyairport.com.au and 6sqft.com. 


The goal of this method is to transfer the main characteristics of the image from WF to VF. The scheme of the 
method is very simple, although the operations have a large computational volume. The method is implemented in five 
stages (Fig. 1): 

1) Selection of color (RGB) channels for WF — WFr, WFg, WFb and for VF — VFr, VFg, VFb. The RGB model 
is used because of its generality and correlation between channels, which is specific to this model. 

2) Calculation of two-dimensional direct Fourier transform (direct Fourier transform — DFT): DFT(WFr), 
DFT(WFg), DFT(WFb), DFT(VFr), DFT(VFg), DFT(VFb). The DFT formula is given below: 

M-\N-1 
XD =F DL Almere jens Ty (l) 
where c = R, G, B — index for red, green and blue color image channels; M, N — row and column numbers of the 
pixel matrix of the transformed image; k, / — spatial frequency arguments; x,(m,n) — pixel value with spatial 
coordinates (m,n) in channel c; X.(k,/) — complex numbers. 

3) This is a key stage. Transplantation of low-frequency part (LFP) is carried out between pairs of WF and VF 
spectra for each of the red, green and blue channels. This means that VF LFP is replaced by the corresponding WF 
LFP. The idea of spectral transplantation is based on the following property of DFT: the general character of the image 
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(i.e., color hue, brightness, contrast) depends on the spatial frequencies contained in LFP (including the constant 
component) of its two-dimensional spectrum. 

Thus, by transplanting WF LFP into VF spectrum, we transfer the main characteristics of the image from WF to 
VF. For this, it is more convenient to use a centered form of a two-dimensional spectrum, where the constant 
component is located in the center of the matrix of spectral coefficients, and the low-frequency components are 
symmetrically arranged around the constant component. In a centered spectrum, LFP is the central part of the DFT 
matrix, and LFP has the size M1 x NI (MI <M, NI <N). If M1 =NI, then the notation for the square matrix LFP can be 
LFP(012..F), where 0 — constant component, F — number of the largest spatial frequency in LFP matrix. 

The size of the LFP for transplantation depends on the size of the transformed image (this size determines the 
spectral resolution) and on the volume of image characteristics that should be borrowed from WF. At this stage of 
research, the size of LFP is determined empirically. For example, the best visual results for 512x512-pixel images 
were obtained using LFP (012345). 

4) Restoring RGB channels for VF using a two-dimensional reverse Fourier transform (RFT). While at this stage, 
the characteristics of WF and VF are mixed. As a result, RGB channels of the VF image are obtained with the main color, 
brightness and contrast characteristics of the WF, as well as with characteristics inherited from the original VF. 

5) Restoring the corrected VF color through merging the RGB channels obtained at the previous stage, cutting out 
virtual objects and building an AR scene by superimposing the cut virtual objects on WF. 

Obviously, if this method is used to process the WF video stream, then there is no need to calculate DFT, and, 
accordingly, LFP for each frame of the real world, since the main characteristics of the image are chan ged only with a 
radical change in the recorded scene. Such changes can be easily detected by jumps in the average pixel value. At these 
moments, it is needed to recalculate the spectral transformation for LFP. 

Since various types of images are optimally described by different types of spectral transformations mentioned 
above, it is reasonable to develop an automatic algorithm for selecting the optimal type of transformation for use in 
spectral transplantation. 

We propose to estimate the difference between the visual perception of VF and WF by the RMS distance 4 between 
the LFP power spectra of their images (for all color channels): 

+S SB. (kD — Pyksl) 


) k=01=0 


Ave (Ree 


L 


(2) 
where Py and Pw — two-dimensional power spectra of VF and WF, respectively. For example, in the case of the 


Fourier transform, the formula for P has the form: 


M-1N=1 km nl. | 
P (k,l) = (m, —j2n(—+—))] . 
(kL) =| XY x, (m,n) exp(-j oT are, 


m=0 n=0 


(3) 


We propose to determine the optimal type of spectral transformation by the proximity of the vectors 4 and the mean 
vector calculated by the criterion of the minimum sum of squares of the distances between the mean vector and the 
vectors A for all transformations under consideration. 

Let 4;(A;R, 4jG, 4;B) be the normalized vector of the distance between spectrum VF and WF LFPs for conversion j. 
Let Aq(AaR, MaG, AaB) be mean vector, and D; — distance between 4; and 4,. Then, the sum of S squared distances 


from the vectors 4; of all transformations under consideration to the mean vector is equal to: 


S=ED) =T TA, -A,), c= RGB. (4) 
P | joe 
Coordinates Agr, 4ac, Jap of the mean vector are calculated as the solution to a system of partial differential 
equations: 
25: c=R,G,B. (5) 
OA,. 


The selection of the optimal type of spectral transformation is determined by the proximity condition: 
min D,. (6) 
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Another obvious criterion for selecting the optimal type of transformation is the length of the vectors 4. However, 
the extremes of such a criterion may be related to the ability or inability of certain transformations to correctly detect 
the difference between certain types of WF and VF. Therefore, we consider the use of the mean vector as a more 
reliable method of selection. 

Similar to the DFT calculation for WF, the optimal type of transformation is selected only once at the beginning of 
spectral transplantation, unless WF is changed dramatically. 

Research Results. The proposed method was tested using the Fourier transform without selecting the optimal type 
of transformation. WF (real airport scene) and VF (virtual airplane model) had a size of 512=12 pixels and 24-bit 
colors. Two different conditions were investigated: 

1) WF — photo of the airport in cloudy weather (Fig. 2 a); 

2) WF — photo of the airport in sunny weather (Fig. 2 b). 

In both cases, VF contained a 3D model of the aircraft shown in Figure 2 c. LP(0), LP(01), LFP(012), LFP(0123), 
LFP(01234), LFP(012345) transplants were tested. Some of the test results are shown in Figure 3. The best visual 
results were obtained using LFP(012345). In Figure 3, the images after spectral transplantation are intentionally shown 


without other VC effects (shadows, lighting, etc.) to demonstrate the pure results of this method. 


LFP(0123) : LFP(012345) 


Fig. 3. AR-scene: a — consisting of WF and VF without LFP transplantation; b — AP-scene after transplantation 
LFP(0123); c — AR-scene after transplantation LFP(012345) 


The upper and lower rows in Figure 3 correspond to the opposite conditions for WF: light and dark WF with 
different shades. Experiments with any other WF will not add significantly new information since they will have 
conditions between those already presented in Figure 3. 

Numerical simulation was carried out to demonstrate the mechanism of spectral transplantation. Figure 4 shows 
Fourier transplantation using a small (8x8) pixel matrix representing one of the color channels WF and VF. Such a 
small size of the matrix enables to clearly illustrate the transplantation procedure. In this example, the WF matrix can 
be associated with an image with a vertical gradient fill, and the VF matrix — with an image with a horizontal gradient 
fill. Another difference between WF and VF is the range of pixel values: 8-15 for WF (“lighter image”, 8 is a constant 
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component) and 0-7 for VF (“darker image”). LFP(01') transplantation is shown, where |' means part of the first spatial 
frequency component (used because of the very low resolution of the 8x8 matrix). The 3D form of the VF matrix 
after transplantation indicates the transfer of properties from WF to VF: the edge of the surface has risen; the first pixel 
has received the value of the constant WF component. This example demonstrates how, as a result of spectral 


transplantation, VF starts to acquire a vertical gradient and a constant component. 


3D view of WF pixel matrix 3D view of VF pixel matrix 3D view of VF matrix with LPF (01’) 
i : 1o aN 20 
0 AZZ IS 8 j <> 15 
La AY LY sr sar a | 6 i i a SFE 
5 ad 4 f 10 LIST aa OS 
tL \ 7 EE SET SEE 7 
0 _2 » NEE T 
0° 0 — 
Two-dimensional WF spectrum Two-dimensional VF spectrum VF spectrum with LFP (01’) 


LFP(1 ae 
LFP(0) LFP(01’) transplantation 


Fig. 4. Numerical simulation of spectral transplantation for 8x8-pixel matrices 


Spectral transplantation provides several options for changing the parameters of this procedure: changing LFP size; 
selecting individual components of the spectrum for transplantation; using different transplant coefficients for various 
components to be transplanted. 

Figure 5 shows the effect of transplantation with different parameters for various types of virtual objects — virtual 
aircraft models that differ in surface texture, markings, and gloss. Figures 5 a and b depict a virtual airplane with 
complex textures, text symbols and reflections of virtual light sources. Figures 5 d, e and f show a virtual plane with 
simple contrasting colors. Parts a and d contain virtual objects without transplantation; b and e contain virtual objects 
after LFP transplantation(0123); c and f contain virtual objects after LFP transplantation(012345). Virtual objects are 
intentionally shown without other VC effects (shadows, lighting, etc.) to demonstrate the pure results of the method. 
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Fig. 5. Scenes with cloudy WF: a, d— AR-scenes consisting of WF and VF without LFP transplantation; 
b, d— AR-scenes composed after LFP transplantation (0123); c, f— AR-scenes composed after LFP transplantation (012345) 


It is important to emphasize that the presented figures illustrate the possibilities of tuning the proposed method, 
and not the final result, since it requires tuning to specific WF. The demonstration of well-executed, but incomplete 
results, as is often practiced in VC works, does not seem correct to us. 

Discussion and Conclusions. The key complicating factor for the described method, presented in Figure 1, is the 
high computational costs. The most promising way to solve this problem is to directly convert WF LFP parameters 
into VF rendering parameters. This eliminates the cumbersome procedures of three DFT and three RFT calculations 
at the second and fourth processing steps, and requires only three WF DFT calculations, once for each section of the 
WF flow without significantly changing the basic characteristics of WF. This approach provides processing the VF 
flow in real time. 

Another problem is selecting the optimal LFP size. As the volume of spatial frequencies used increases, they begin 
to hold information about the WF contents. Therefore, limiting the size of LFP is needed to eliminate the effect of a 
hybrid image [11]. The complexity of the optimal selection is conditioned by its association with both the LFP size 
and the nature of the image. Recent advances in deep learning suggest that a new approach related to visual coherence 
through spectral transplantation could be the use of generative adversarial networks (GAN) to transmit realistic lighting 
information from the source image to the target image in the same way that GAN do to transmit image style. In 
particular, it would be interesting to compare the performance of GAN in the case of data sets consisting of either RGB 
images or images represented in the frequency domain using DFT. We believe that the latter approach will help to 
select the optimal LFP. GAN are already widely used in VC study [2] as are neural networks in general [12]. 

In further research related to the topic of this paper, the following issues will be considered: 

— automatic determination of the optimal LFP size for transplantation with a given volume of characteristics 
borrowed from WF; 

— automatic detection of the exact moments when it is required to calculate new WP LFP for transplantation when 
processing WF video in real time (as mentioned above, this must be done if the basic characteristics of WF are radically 
changed); 

— using the same approach in reverse (from virtual to real) to apply virtual lighting to real scenes (how virtual 
lighting affects the environment). 

As a fully automatic process without measuring illumination, the proposed spectral transplantation method solves 
a number of complex VC problems. Let us say, how to best align the color, brightness, and contrast characteristics 
between real and virtual components in AR scenes. All these tasks are solved through one simple procedure without 
modeling lighting conditions, AR-scene geometry or BRDF, which eliminates the inevitable modeling errors. The 


proposed method can be a valuable addition to other VC tools. 
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Abstract 

Introduction. It is obvious that in the near future, the issues of equipping moving robotic systems with autonomous 
control elements will remain relevant. This requires the development of models of group pursuit. Note that optimization 
in pursuit tasks is reduced to the construction of optimal trajectories (shortest trajectories, trajectories with differential 
constraints, fuel consumption indicators). At the same time, the aspects of automated distribution by goals in group pursuit 
were not considered. To fill this gap, the presented piece of research has been carried out. Its result should be the 
construction of a model of automated distribution of pursuers by goals in group pursuit. 

Materials and Methods. A matrix was formed to study the multiple goal group pursuit. The control parameters for the 
movement of the pursuers were modified according to the minimum curvature of the trajectory. The methods of pursuit 
and approach were considered in detail. The possibilities of modifying the method of parallel approach were shown. 
Matrix simulation was used to build a scheme of multiple goal group pursuit. The listed processes were illustrated by 
functions in the given coordinate systems and animation. Block diagrams of the phase coordinates of the pursuer at the 
next step, the time and distance of the pursuer reaching the goal were constructed as a base of functions. In some cases, 
the location of targets and pursuers was defined as points on the circle of Apollonius. The matrix was formed by samples 
corresponding to the distribution of pursuers by goals. 

Results. Nine variants of the pursuit, parallel, proportional and three-point approach on the plane and in space were 
considered. The maximum value of the goal achievement time was calculated. There were cases when the speed vector 
of the pursuer was directed arbitrarily and to a point on the Apollonius circle. It was noted that the three-point approach 
method was convenient if the target was moving along a ballistic trajectory. To modify the method of parallel approach, 
a network of parallel lines was built on the plane. Here, the length of the arc of the line (which can be of any shape) and 
the array of reference points of the target trajectory were taken into account. An equation was compiled and solved with 
these elements. On an array of samples with corresponding time values, the minimum time was found, i.e., the optimal 
time for simultaneous group achievement of multiple goals was determined. For unified access to the library, the control 
vector was expressed through a one-parameter family of parallel planes. A library of calculations of control vectors was 
formed. An example of applying matrix simulation to group pursuit was shown. A scheme of group pursuit of multiple 
goals was presented. For two goals and three pursuers, six samples corresponding to the distribution of pursuers by goals 
were considered. The data was presented in the form of a matrix. Based on the research results, the computer program 
was created and registered — “Parallel Approach on Plane of Group of Pursuers with Simultaneous Achievement of the 
Goal”. 

Discussions and Conclusion. The methods of using matrices in modeling group pursuit were investigated. The possibility 
of modifying the method of parallel approach was shown. Matrix simulation of group pursuit enabled to build its scheme 
for a set of purposes. The matrix of the distribution of pursuers by goals would be generated at each moment of time. 
Methods of forming matrices of the distribution of pursuers and targets are of interest in the design of virtual reality 
systems, for tasks with simulating the process of group pursuit, escape, evasion. The dynamic programming method opens 
up the possibility of automating the distribution with optimization according to the specified parameters under the 
formation of the matrix of the distribution of pursuers by goals. 
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AHHOTauHA 

Beedenue. OuesugqHo, 4uTo B OnMxKalilIee BpeMA COXpaHAT aKTYaJIbHOCTb BOMPOCbI OCHAaLeHHA ABYKYLIMXCcA 
PoOoTOTeXHHYeCKHX KOMIIICKCOB 3JIEMCHTaMH ABTOHOMHOLO yupaBJIeHHA. ITO TpeOyeT pa3sBHTHA MOJeNe TpymmoBoro 
lipecueqoBanua. OTMeTHM, YTO ONTHMUB3AaIlHA B 3aa4yax TpeceqOBAaHHA CBOAHMTCA K MOCTPOCHHIO OMTHMAJIbHBIX 
TpaeKTopuit (KpaTuaiwiMe TpaeKTOpuH, TpaeKTopuu c AuddPepenuHabHEIMH OTpaHH4eHHAMH, NWOKa3aTeIH pacxoya 
TormBa). [Ipu 9TOM He paccMaTpHBaIOTCA ACIICKTLI ABTOMATH3HPOBAHHOLO paciipeseeHHA M0 WesA1M Ip rpyiioBoM 
mpecneqoBanun. J|ni BOCHOMHeEHHA ITOTO Tpobea BEIMIONHeHa MpescTaBeHHat HaydHaa padota. Ee pesynbTaToM 
JOJDKHO CTaTb MOCTpoeHHe MOJeIM ABTOMATH3HpOBaHHOrO paciipeyeueHuA MpecteqoBaTeei M0 WesAM B IpyilloBoM 
TipecleqOBaHHun. 

Mamepuaanoi u memoooi. [na w3ysenus rpymmoBoro MpecweqoBaHuA MHOXKeCTBa ese ctopMHMpoBaHa MaTpHa. 
YmpaBiAloWve MWapaMeTph! ABUWXKCHHA MpecweqoBaTeseH MOAMPUUMpOBaHbI HO MMHUMaJIbBHOM KPHBH3He TpacKTOpHH. 
J\eTanbHo paccMOTpeHbI MeTOAI MOroHH U cOmMKeHUA. [loKa3aHbI BOSMOKHOCTH MOAMMUKALMM MeTOsa MapaseubHOrO 
cOmmxkeHu1. MatpwuHoe MoyesMpoBaHve 3ajelcTBOBaIM JJId MOCTpOeHHA CXeMbI TpyiMoBoro MpecweqOBaHHA 
MHOxecTBa Hemel. [lepeuucueHHble MpOLecchl MpOMWIOCTpHpoBaHbl (yHKUMAMM B 3aJjJaHHbIX CHCTeMaX KOOPAHHAT U 
aHuMaunei. Kak 6a3a (yHKUMi MOCTpOeHbI OIOK-CXeMBI (Pa30BbIX KOOPAHHAT MpecieqoBaTelA Ha CHeYIOWeM Ware, 
BPeMeHH HM paccTOAHHA JOCTIMKeEHHA MIpecseqoBaTesem Wer. B paye ciryyaeB pacnouoxKeHHe Lelel u mpecieqoBatesen 
ompeyeseHo Kak TOUKH Ha OKpyxkHocTH AnosioHHa. Matpuua ccbopmupoBaHa 0 BbIOOpKaM, COOTBeTCTBYIOUIMM 
paciipeyleseHuro pecseqoBaTesen M0 IeuaM. 

Pezyismamoi ucciedoeanua. PaccMoTpeHbl JeBATb BapHaHTOB MOFOHH, MapasWIeNbHOrO, MpOMOpyMoHabHorO 
TPeXTOYCYHOTO COMM KEHHA Ha MIOCKOCTH H B MpoctpaHcTBe. PaccuvTaHo MaKCHMAaJIbHOe 3HayeHHe BPeMeH JOCTWKeHHH 
lee. OTMe4eHBI CyyaH, KOra BEKTOP CKOPOCTH MpecieqoBaTeIA HalpaBJIeH MPOU3BOJIbHO H B TOUKY Ha OKpy?KHOCTH 
AnlosIoHHa. OTMeYeHO, YTO TpexTOYeYHBI MeTOA cOMMKeEHHA yHOOeH, eCIM Web TBMKeETCA MO OaMcTHYecKOH 
Tpaextopuu. Jia MoAMuKayHH MeTOAa NapaliesbHoro COMKeHHA Ha WIOCKOCTH CTPOHTCA CeTb MapasIebHbIx 
muni. [pu 3TOM y4TeHbI WHHa Ayu MH (KOTOpad MOET ObITb MIPOH3BOJIbHOM (POpMBbI) H MaCCHB OMOPHBIX TOUCK 
TpaekTopuu wel. C aHHbIMH 93JIEMeHTaMHM cCOcCTaBIeHO HM pelleHoO ypaBHeHue. Ha maccuBe BbIOOopoK c 
COOTBETCTBCHHBIMH 3HaYCHHAMH BPeMeH HaliJeHO MMHMMaJIbHOe BPeMA, TO eCTb OMpeeeHO OMTHMasIbHOe BpeMA 
OJHOBPeMeHHOTO TpyMMoBoro OcTwWKeHHA MHOxecTBa Lene. Jai yHudpuyupoBaHHoro oOpaujenua K OnOMOTeKe 
BbIPAKEH YIPaBIIAIONIMM BEKTOP Yepe3 OAHOMApaMeTPHYECKOEe CeMEMCTBO MapaWIeIbHBIX WWIOcKocTeli. CpopMupoBana 
OuOMOTeKa pacueTOB yiIpaBJIAIOWMx BeKTopoB. Hloka3aH UpHMep UpHMeHeHHA MaTpH4Horo MOJeMpoBaHHA K 
rpylmoBomy lpecsieqoBannio. Ipeyctapsena cxeMa rpymoBoro NpeceqoBaHuaA MHOxKecTBa Leen. Ja WByX Weevi u 
Tpex IpecuieqoBaTesei paccMOTPeHbI LWeCTb BLIOOPOK, COOTBETCTBYIOINMX paclipeseeHHIo IpecieqoBaTesen M0 WesAM. 
JlanHble npeycTaBeHbI B Bue MaTpuubi. Ilo uToraM Hay4HbIx U3bICKAaHH Co3qaHa UW 3aperMcTpupoBaHa UporpamMa 
aia IBM «Moyen mapanenbHoro cOmMWKeHHA Ha MWIOCKOCTH IrpylIlibl WpecweqoBaTeseli Cc OHOBPeMeHHbIM 
TOCTWKEHHEM Les)». 

O6cyscdenue u 3aKkniouenue. VUiccneqoBaHbl MeTOAbI MCHONb30BaHHA MaTpHI, pH MOJeMpoBaHuu rpynnmoBoro 


TipeciI¢qOBaHHaA. Tloka3aHa BO3MO2KHOCTb Moqu@uKalMn MeTOJa = TlapaJIeIbHOTO cOIMKeHUA. Marpyu4Hoe 
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MOJCJIMpOBaHHe TpymIOBoro WpecsIeqOBaHHA WO3BOJIACT BbICTPOHTb CrO CXCMYy JIA MHO2KCCTBA wWeseu. Marpnia 
paciipeqesIeHuA lipecsIeqOBaTesen TIO WeJIAM OyqeT TeCHepHpOBaTbCsA B KaK bIM MOMCHT BpeMeHH. MeroybI cbopMupoBaHHua 
MaTpHI paciipeqeseHua lipecuieqoBaTesieit Mi Ween TIPCACTaBILAFOT HHTepec Ip WpOeCKTHpOBaHHu CHCTCM BUpTyasIbHoHi 
peasIbHOCTH, It 3afad C MOCIMpOBaHHeM Iipowecca rpylHOBoro IipecsiIeqoBaHHA, yOeraHua, YKJIOHCHHA. Metoy 
AWMHAMHYeCKOrO HpOrpaMMHpoBaHHa ITIpu cbopMupoBaHHH MaTpHUbI paciipewqesienuA lipeculeqoBaTesei TO WeJIAM 
OTKPBIBaeT BO3MO?2KHOCTb aBTOMAaTH3alHn paciipeesIeHHA C ONTHMH3alMe 10 3a{aHHbIM TlapaMeTpaM. 


KoroueBbie  cJIOBa: ayTOpHTM IpyHMOBOrTO HpecsiIeqOBaHHA, ONTHMH3alvaA B 3ayadax = lipeciueqOBaHHA, 
aBTOMAaTH3HpOBaHHoe paciipewesenne TlO IeJ1AM, MaTpHia AOCTH2KCHUA TipeciIeqOBaTeAMH lWeseH, 
aBTOMAaTH3HpOBaHHoe NpHHATHe pemiennii, aBTOHOMHOE yiipaBJIeHHe, MapasIesIbHOe cOmmKeHHe, TIPONOPUHOHAIBHOe 
cOmmKeHHe, TpexToueyHbli MeTO cCOIMKeEHHA, OuOmMoTeKa pacudeToOB yiipaBJIAIOWIHX BEKTOPOB 


BaarogqapHoctTn: aBrop BbIpaxkKaeT MpH3HaTebHOCTS ZupeKtopy MHctuTyta MaTeMaTHKH Hu UHdopMaTuKu bypatcKoro 
rocyfapcTBeHHoro yHuBepcnteta uM. JJ. ban3apopa AnToHoBoi Jlapuce BacusbeBHe 3a MOMOLIb, OKa3aHHy!0 B padote 
Hay, cTaTbeii. 
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Introduction. The pursuit algorithms are studied from the point of view of their classical and optimal implementation. 
Their role in differential pursuit games is investigated. The applied sphere of ready-made solutions is very wide, because 
the results of such scientific research are applicable in various information technologies and systems, in particular, in 
search engines. Undoubtedly, the issues of equipping moving robotic complexes with autonomous control elements will 
be of current concern for a long time, which also requires high-quality implementation of the algorithms under 
consideration. 

In [1-4], the coordinated behavior of a group of pursuers and targets was investigated. For general theoretical and 
practical issues in the problems of persecution, works [5—9] were considered. The guidance of the pursuer to the target 
was analyzed considering the information provided in [10-13]. 

With all the theoretical and practical interest in this topic, optimization in pursuit problems was limited to the 
construction of optimal trajectories. Specifically, the shortest trajectories, trajectories with differential constraints, fuel 
consumption indicators were proposed. But the aspects of automated distribution by goals in group pursuit were not 
considered. To fill this gap, this scientific work was carried out. Its basic result was the construction of a model of 
automated distribution of pursuers by goals in group pursuit. The formation of a matrix of achieving goals by pursuers 
was shown. When assigning goals to the pursuers, all possible combinations of achieving goals were sorted out, and a 
combination of the minimum value of the criterion from the generated set with the maximum value was selected. 

Optimization of the multiple goal group pursuit is a promising direction for the development of such a discipline as 
optimal motion control in tasks related to automated decision-making and autonomous management. 

Materials and Methods. In the model of group pursuit described in the paper, targets move along predefined 
trajectories. However, this predestination does not matter in principle. The pursuers are distributed by the targets 
automatically, based on the minimax solution of the goal function. Then the control parameters of the pursuers' movement 
are modified. In this paper, this is the parameter of the minimum curvature of the trajectory. This approach allows for 
simultaneous achievement of goals. 

Consider a group pursuit of a set of goals: N pursuers catch up with M goals. We form a matrix of the distribution of 
pursuers by goals: 

YW, where i=1.N,j=1.M. 


Each cell ‘Yj; contains information about the phase coordinates of the i-th pursuer and the j-th target. Matrix Vj 
contains information about the method by which the i-th pursuer goes after the j-th goal. 

The data stored in the cells of the matrix determines the access to the library of calculations of the control vectors of the pursuer. 

In each cell of matrix Yj, the predicted time for the i-th pursuer to reach the j-th goal can be calculated: fj. 

Research Results 


In each received sample A, ={v hs it is required to find the maximum value of achievement times 


Adie 77 ~~ tnk Ink 


i Max{t,| ,e.g., from {t,,, (Table 1). 
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Table | 
Samples corresponding to the distribution of pursuers by goals 
Goals 

g 2 1 2 1 2 1 2 1 2 1 
zB 1 x x x 
5 
a 2 x x x x 

3 x x x x 
Samples Al A2 A3 A4 As A6 
It is necessary to form matrices ‘Yj, where i= 1...3, j=l...2 according to the possible samples 


A,k=1...6. Then, after the conversion, maximum value t;=Max {t,} is found. The calculation made it possible to establish 
that pursuer P;, demonstrated the greatest time of achievement, catching up with goal 7; from sample Az. 
Thus, consider sample Ax. You can increase to the value of parameter ¢, all values tj, depending on the velocity vectors 


of the pursuers and goals, as well as their permissible angular velocities. This determines maximum value t. 


Having received an array of samples {Ax} with corresponding time values {tj}, it is necessary to find minimum time 


tmin=Min { t,}. This is how the optimal time for simultaneous group achievement of multiple goals is determined. 
Algorithms for calculating the next step of the pursuer and estimating the time when the pursuer reaches the 
goal. Figure 1 shows the algorithm of the function for calculating the next step and the speed vector of the pursuer. 


Accessing the library of control vector 
Start age 
calculations u 


Input of current coordinates of pursuer P , vector of the current speed of 


pursuer V,, , permissible angular rotation of pursuer wp, control vector 7, 


Pe 
discrete time interval Ar. 


Translation to local coordinate system. Abscissa is vector # . Start is at the 
location of pursuer P. Speed vector of pursuer V, . In the local coordinate system, 


angle a between velocity direction V, and control vector # (abscissa of the local 


coordinate system) is calculated. 


‘OS ( Oey ) 
IM( Gey) 


Translation of values P,,,,,V,, into the Translation of values P.,,,V,,.,, into the 


new?" Ihew new ? 


traditional coordinate system. traditional coordinate system. 


Fig. 1. Flowchart for calculating phase coordinates of the pursuer in the next step 
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Figure 2 shows an algorithm for calculating the time and distance of the pursuer reaching the goal. Variable e€ is the 


threshold value of the distance from the pursuer to the target, at which the goal is considered to have been achieved. 


Calculating the distance between pursuer Entering the initial phase coordinates 
and target AS, i = 0. of pursuer and target. 


Reference to function of 
calculating phase 


Calculating the distance between pursuer and coordinates of pursuer. 


target AS, i=i+1. 


Entering current phase 
coordinates of target. 


Fig. 2. Flowchart of the function for calculating the time and distance of reaching the goal by the pursuer 


If the target moves along a predetermined trajectory, then the algorithm shown in Figure 2 can give an estimate of 
time f; for the i-th pursuer to reach the j-th goal. In this case, the output parameter of the function can be the number of 
iterations of the pursuit process Nj; Number of iterations Ni;— output parameter of the function for calculating the time 
and distance of reaching the goal by the pursuer. 

If the goal replies to avoid achievement, you should evaluate the time differently. It is necessary to build predicted 
trajectories as composite segments of straight lines, arcs of circles, square and cubic parabolas and other known lines. 
This will make it possible not to solve boundary value problems in the calculation cycle. 

Formation of a library of control vector calculations. The distribution matrix ‘i, where i= 1...N, j=1...M 
pursuers by goals is built on each discrete time interval. In each cell of matrix Yj; , information about the method 
of persecution is stored. It is based on the reference to the library of functions for calculating control vectors u 
(Fig. 3-11). 


Fig. 3. Pursuit method on the plane and in space: 4; = 


. Here, T, — target position point, 


P — point of the pursuer's position [14] 
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|K-P| 


Fig. 4. Method of parallel approach on the plane: u = T — target position point, P — point of the pursuer's position, 


K — point on the Apollonius circle. It is uniquely determined by points P, T and target velocity vector V, [15] 


|K-P 


on the Apollonius circle. The Apollonius circle lies in plane X, 


Fig. 5. Method of parallel approach in space: u = , T — target position point, P — point of the pursuer's position, K — point 


formed by points P, T and target velocity vector Vr{16]. The case is shown when the velocity vector of the pursuer is directed 


arbitrarily. As time passes, the velocity vector of the pursuer is directed to a point on the circle of Apollonius. 


z=fQy). 


Fig. 6. Pursuit method on the plane: u, = 1 , where P,, — result of the intersection of surface z=f (x,y), planes Ae and 


Fi, 
IF. 


. 


sphere S; centered at point Pi. Paguyc \v,,| -At . Radius ra — orthogonal projection of point P; onto XY plane. 


For unified access to the library, it is necessary to express the control vector [17]. 
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i+] i 


Fig. 7. Method of parallel approach on the plane: u, = . Here, Pi+1 — result of the intersection of surface 


z=f (xy), plane PoP T,,, and sphere Sicentered at point Pi. Radius \v,,| “At. 


i+” itl 


Point ae — orthogonal projection of point Pi+1 onto XY plane. For unified access to the library, it is necessary to express the 


control vector. ®;— one-parameter family of parallel planes [18] 


d0_ dy FAM le 7 
Fig. 8. Proportional approach method: — = k -— , Ag = arccos ; 
dt dt 2 ri T| : T.., 
Beha aa tal 
; ; , V,-At-cos| k-arccos 

inf +nal In -T, IE) Ke _ pp 
A0 =k-arccos| — ; —— a ae =— : 
2 ARI WE Fal - Tra i 


nl + 
V, -At-sin| k -arccos 5 


T|- 
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Fig. 9. Three-point approach method: 


The method is convenient if the target is moving along a ballistic trajectory. 
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fia Cs) 


f(s) 


Fig. 10. Modification of the parallel approach method on the plane. A network of parallel lines is being built 
Tig (s) =f, (s) +T,,,—T, , where s — arc length of the line, T;— array of reference points of the target trajectory. Solving the 


equation (ha (s) = py = (V, -Aty with respect to parameter s provides finding value s", which corresponds P., = f,,, (s° ). 


-—P 
u, = [P.,-P] Family fi(s) can have lines of any configuration [19]. 


Fig. 11. Modification of the parallel approach method on the plane. Network /i(s) is built, where s — arc length of the line, 
T; — array of reference points of the target trajectory. The condition is fulfilled that the end of line fi(s) passes through point T; and 


point P; is an incident on line fi(s), i.e., it is used as a pattern. Solving the equation (fa (s) —P y = (V, -Aty with respect to 


& 4 tan F 
parameter s provides finding value s*, which corresponds P,,, = f,,, (s ) JU, = iP. —P] : 


i+] 


Family fi(s) can have lines of any configuration [20]. 


Thus, the library of calculations of control vectors contains methods of pursuit on the plane, in space, and on the 
surface. Parallel approach methods are calculated on the plane, in space, and on the surface. Proportional approximation 
methods are calculated on the plane and in space. Three-point methods are calculated on the plane and in space. Modified 
pursuit methods are calculated on the plane and in space, when the permissible curvature of the trajectories is used to 


control the pursuer. Modified methods of parallel approach are calculated on the plane and in space. 
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Modification of the methods of parallel approach and pursuit provide building a network of predicted trajectories that 
allow for various boundary conditions. This is illustrated in Figures 3-11. But not all methods of calculating control 
vectors are presented in them. It is assumed that this is an open, replenished library of functions. 


Case of applying matrix modeling to group pursuit. Consider a case of group pursuit (Fig. 12). 
I 


Fig. 12. Scheme of multiple goal group pursuit 


Here, all pursuers achieve the goal using a modified method of parallel approach, which corresponds to Figure 10. In 
the pursuit model in Figure 10, the curvature of the trajectory should not be greater than a certain value. Therefore, the 
initial radius of curvature of the trajectory increases for pursuers Pz and P3as shown in Figure 12. 

Sample Ax, has been formed, in which pursuer P; catches up with 7;, Then there is a primary evaluation of the time of 
reaching tj. To estimate time fy, the following are calculated: 

— length of the rectilinear section to the target, 

— length of the arc of the mating circle of the permissible radius. 

Then maximum value 4=Max{tj} is selected. An increase in time fj to % occurs in this model due to an increase in 
radius of the mating circle from value 7; to 7; +67;,in pursuer P;. 


Figure 13 is supplemented with an animated image showing the process of multiple goal group pursuit [21]. 


Trajectories of pursuers Positions of pursuers 
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Trajectories of pursuers Positions of targets and pursuers coincide 


-30 = = & 
-30 30 


Networks of predicted trajectories Positions of targets 
b) 
Fig. 13. Schemes of group pursuit phases: a — initial phase; b — final phase 

Based on the results of the study, a computer program, which implements an algorithm for group pursuit of several 
goals was created and registered [22]. This software solution is called “Parallel Approach on Plane of Group of Pursuers 
with Simultaneous Achievement of the Goal”. 

Discussion and Conclusions. The methods of pursuit, parallel, proportional and three-point approach were described 
and visualized as functions on the plane, on the surface, and in space. In addition, the possibilities of modifying the 
method of parallel approach on the plane were shown. With the application of matrix modeling to group pursuit, a scheme 
of multiple goal group pursuit was built. The initial and final phases of this process were shown separately. The calculation 
of the achievement time made it possible to identify the pursuer who needed the most time to reach the goal from the 
sample under consideration. 

Thus, it is assumed that the matrix of the distribution of pursuers by goals is generated at each moment of time. Goals 
and pursuers may disappear, new ones may appear. This matrix can also be used by the party representing the targets who 
evades prosecution. The results of the scientific research described in the article allow us to form the principles of 
automated distribution of pursuers by goals based on the selected target function. Algorithms for modifying the 
trajectories of pursuers to achieve goals simultaneously or according to a set schedule were proposed. The issues of 
forming a library of pursuit methods were also considered. The method of forming a matrix of the distribution of pursuers 
by goals can be in demand when designing virtual reality systems for game tasks in which the process of group pursuit, 


escape and evasion is simulated. 
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Abstract 

Introduction. The challenges of placing virtual objects in a real-world environment limit the potential of augmented 
reality (AR) technology. This situation identifies a gap in scientific knowledge that requires additional research. 
Therefore, the main task of this study was to develop a method for optimal placement of virtual objects, in which the 
objective function of comfort was minimized. This approach is aimed at improving AR systems and developing the 
corresponding theory. 

Materials and Methods. The conducted research was based on the analysis of the placement of virtual objects in AR/VR 
applications with particular emphasis on optimization. The concept of comfort of placement was proposed, taking into 
account the size of the object and the distance to the boundaries of free space in X, Y, Z coordinates. 

Results. As part of the study, formulas were obtained for the optimal placement of objects with an arbitrary comfort 
function. The basic criterion was to minimize the difference between comfort levels from different sides of the object. It 
was found that a successful placement of objects required taking into account their size and comfort zones, as well as 
solving a system of n linear equations. 

Discussion and Conclusion. The results obtained make an important contribution to the study of the problem of placing 
virtual objects in AR/VR/MR. They open up new opportunities for improving user interaction and conducting further 
research in the field of spatial computing. Possible directions for further development are dynamic adjustments and 
integration of the results into various XR scenarios. 


Keywords: augmented reality, virtual objects, physical space, optimal placement, mathematical model, equations 
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Hayunaa cmamoa 


Pa3mMeleHHe HECKOJIBKHX BUPTYaJIbHbIX OOLCKTOB B PU3H4eCKOM 
IIpOCTpaHcTBe B NPH.IOKCHHAX JONOJHEHHOM peasIbHOCTH 
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Mockosckuit NOIMTexXHuYeCKHH yHuBepcurter, r. Mocxsa, Poccuiickaa Dexepayua 
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AHHOTauHA 

Beedenue. \[poOnempl, cBa3aHHble Cc pa3MelleHHeM BHpTyaIbHbIX OOBEKTOB B PpeasIbHOM cpeye, CyLeCTBeHHO 
OrpaHH4HBalOT BO3MOXKHOCTH TeXHOJOrHH AOMONHeEHHOM peambHocTH (AR). Takaa cuTyallua BbIABJIAeT TIpOOeN B 
Hay4HBIX 3HaHUAX, TpeOyIOMMi JONOJHUTeIbHOrO UcceqoBaHHA. IlosToMy OCHOBHOH 3ayauel JaHHoro HcceqoBaHUA 
ABHIACbh pa3paOOTKa MeTOa OMNTHMaJIbHOrO pa3MeCHHA BUPTyaIbHbIX OOBCKTOB, Tp KOTOPOM MpPOHCcxXOAHT 
MHHUMH3alHA WeseBon PyHKUMH KOM@opTHocTH. Tako WoAxXoy HalpaBleH Ha yCOBepuIeHCTBOBaHHe cuctem AR u 
pa3BHTHe COOTBeETCTBYIOMeH TeOpuHu. 

Mamepuaaoi u memooui. Uipopeyennoe uccreqoBaHve OCHOBbIBaCTCA Ha aHaJIN3e Pa3sMeLLCHHA BUPTYAIbHbIX OOCKTOB B 
AR/VR mpysioxKeHHAx C OCOOBIM aKICHTOM Ha ONTHMH3alHI0. bpi0 mpes10oxKeHO MOHATHE KOM@OPTHOCTH pa3MelleHHA, 
YUUTbIBAaIOLee pasMepbI OObEKTa H paccTOAHHA JO TpaHHL CBOOOAHOrO IpocTpaHcTBa 0 KoopAuHaTam X, Y, Z. 
Pezyismamot ucciedoeanua. B pamkax vccieqoBaHHaA OLIN MoyYeHbI POpMyBI [IA ONTHMAabHOTO pa3MeLleHHA 
OOBEKTOB C IPOM3BOJIbHOM PYHKUMel KOMPoOpTHOCTH. OCHOBHBIM KPHTepHeM ABIIACTCA MMHUMV3AUHA pasHULbl MEK LY 
YPOBHAMM KOM@OPTHOCTH C pa3HBIX CTOPOH OObeKTa. bp0 BbIABIICHO, UTO yCIelWHOe pa3sMeleHHe OOBeEKTOB TpeOyeT 
yueTa HX pa3MepoB H 30H KOMOPTHOCTH, a TakoKe PeLICHHA CHCTEMBI V3 N JIMHeMHBbIX ypaBHeHHH. 

O6écyordenue u 3akmou4enue. TonydeHHble pesyibTaTbI MpelcTaBaioT coOol BaxKHbIM BKIay B HCcueqoBaHHe 
TIpOOJeMBI pa3sMeLleHHA BUPTYAIbHBIX OObeEKTOB B AR/VR/MR. Onn OTKpbIBaIOT HOBbIC BO3MOXKHOCTH JIA yIyYWeHHA 
B3aHMOJIeHCTBUA C TOJb30BaTeIAMH UM MpOBeeCHHA MabHeMUIMX UcceqOBAaHHHM B OONaCTH MpocTpaHcTBeHHbIxX 
BBIYHCICHHM. Bo3MO2KHbIMH HallpaBICHHAMNH JIA WaIbHeMWero pa3sBUTHA ABJIAIOTCA WHHAMMYECKHe KOPpeKTHPOBKH HU 
MHTerpalla MOUYYCHHBIX Pe3YJIbTATOB B pa3sM4HbIe XR-cyenapun. 


KoroueBbie CJ10Ba: AWOHOJIHCHHAA PeasJIbHOCTb, BHUPTYaJIBHbIe OObEKTHI, dbu3sn4eckoe TIPOCTpaHCTBO, pallHOHaJIbBHOe 
pa3MelieHne, MaTeMaTHYeCKadA MOJCJIb, YDAaBHCHHA 


BuaroqapHoctH: JaHHoe UccreqoBaHHe OcyLecTBIeHO Oarofapa (PuHAaHCOBOH NoAWepxKKe POOH B pamKax Hay4Horo 
mpoexta Ne 21—510—07004. Bripaxxkaem npu3HaTeIbHOCTb KOJWIeraM, YYaCTBYIOWIMM B JaHHOM rpaHTe, 3a UX IeHHBIi 
BKJIaJ{ B padoty. Kpome Toro, xoTHM HoOaroapuTb peqakKWUHOHHY!O KOMaHJly %KypHasla U pelleH3eHTa 3a KOMIIETCHTHYIO 
9KCHeEpTH3y HW WeHHbIe peKOMeHallHH M0 yJIy4WIeHHI0 CTaTbH. 
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Introduction. Rapid development of augmented reality (AR) technology opens up new opportunities in various 
fields — from entertainment to education and industrial applications [1, 2]. However, despite considerable achievements, 
there are numerous problems that need to be solved, specifically, in the context of placing several virtual objects in a real- 
world environment. One of such problems is the optimal placement of virtual objects in augmented reality applications 
to provide for optimal and comfortable user experience [3-5]. 

This problem arose in connection with the need for the device to understand physical space. For effective placement 
of virtual objects in the real world, the application should be able to correctly interpret the material environment in which 
the user is located applying sensors and cameras of mobile devices [6]. 

This article presents a new approach to determining the optimal placement of virtual objects in physical space. This 
problem has some similarities with another close topic of generative contextual scene augmentation (CSA), where the 
key goal is to create a harmonious and convenient interaction between virtual and physical objects [7, 8]. However, the 
approach proposed by the authors differs from the one mentioned, since it focuses on determining the optimal distance 
between objects using a monotone comfort function, while CSA takes into account the semantics of the scene, the context 
and the meaning of virtual objects. 

Existing approaches to solving the problem of rational location of virtual objects are usually limited by assumptions 
about the form of the comfort function, and they do not always guarantee the optimality of the solution. This article 
proposes a new approach that is significantly different from those currently used. This makes it more flexible to find a 
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rational arrangement of virtual objects and provides a universal solution for various scenarios and conditions. The term 
“optimal placement” is used for the following reasons: firstly, the final choice in the proposed recommendations is still 
at the discretion of the user; secondly, despite the fact that in the context of this work, the task of minimizing the objective 
function is solved, it contains elements of fuzzy sets. 

Within the framework of this work, the concept of optimal placement of a set of virtual objects in a one-dimensional 
physical space is introduced. A model has been developed that allows solving the problem of optimal placement of a set 
of objects and results, regardless of the choice of the type of comfort function, in a system of linear equations for 
determining the optimal distances between objects. As an example, the solution to the problem for the case with two 
virtual objects is given. 

The objective of this article is to propose and demonstrate a new approach to determining the optimal placement of 
virtual objects in physical space, to establish its applicability and efficiency. Mathematical formulations and methods for 
solving the optimal placement problem are presented, as well as examples of practical application of the results obtained. 
This will show new opportunities for improving the interaction between virtual and physical objects, as well as contribute 
to the development of the theory and practice of augmented reality. 

Thus, this article is aimed at deepening the understanding of the problem of optimal placement of virtual objects in 
physical space and offering a new approach to its solution. The results of the study can be used to create more efficient 
and user-friendly virtual reality systems, as well as for further development of theory in this area. 

Materials and Methods. Placing virtual objects in a real physical space is a task that arises in almost every AR/VR 
application. For all its simplicity, it can create challenges in case of insufficient attention to the issue of optimizing the 
placement of such objects, up to the complete refusal of a large number of potential users to work with the mentioned 
applications. The optimization problem is most acute when it is required to place several virtual objects at once in a given 
physical space. At the same time, even in the case of placing only one object, only recently the concept of “comfort” of 
its placement was formulated and a corresponding model was proposed [9], consisting of the following. 

An object embedded in a three-dimensional physical area is presented as a rectangular parallelepiped with 
characteristic dimensions: / — length; d — width; h— height. At the same time, for each of X, Z, Y coordinates, the 
following concept of placement comfort is introduced. It is clear that the size of the free space should not be less than the 
size of the object, but, in addition, for each coordinate, the concept of comfortable distances from the object to the 
boundary of free space is introduced. For example, for X coordinate, we introduce the concept of comfortable distance on 
the left — D_ and right — D, and, respectively, left and right comfort —— K_ and K, . Denote the distances to the left 


and right of the object to the boundary of free space X_ and X,. We assume that K_ =1,if X_2D_ and decreases to 


0 when approaching zero. For example, for simplicity, let us take linear dependences K_ (Xx a) D_) and K, (xX, / D,) : 


xX 
—, xX <D, 

K_=3D_ . (1) 
1,X >D. 


Dependence K, (xX, / D,) is similar (1). In exactly the same way, we introduce the concept of comfort on one side 


and on the other for Z and Y coordinates. 
If the size of the free space horizontally is L 2 D_+/+D, , then the problem of comfortable placement (K_ = K, =1) 


does not arise; and all problems appear when /<L<D_+/+D,. In this case, the concept of comfort of the object 


placement is introduced, when comfort on the one hand is not obtained at the expense of comfort on the other hand. We 
introduce the target function of comfort: 


K, =(K_-K,) . (2) 
By optimal placement, we will understand such placement, in which minimum K, is achieved. Obviously, this 
happens if K_=K,. 
As it was shown in [9], if dependences K_ (X 7 / D_) and K, (x im / D,) are linear, the minimum of the target function 
(2) corresponding to the optimal placement of the object is attained at the following values X_ and X,: 
D 


BE i 
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X,=L-X_-1=(L Na (3) 


Formulas (3) are simple and convenient for optimal placement of a single virtual object. As noted in [9], their 


disadvantage is that they were derived for the case of linear (1) comfort functions K_ (X_ Fs D.) and K, (xX, / D,) . We 
show that they are always valid if functions K_ (X_ /D_) and K, (X, /D,) are the same function k(x), satisfying the 


condition that it increases monotonically for 0< x <1, and is equal to 1 at x >1. The type of such function 


K_ (X_/D_) = k(x) is shown in Figure 1. 


1.0 


0.5 


Comfort criterion K- 


0.5 1.0 1.5 


Distance normalized to the comfort zoneX- / D- 


j=) 


1-x°,0<x<l 


Fig. 1. Comfort function dependence k(x) = 
l,x>l 


Figures 2 and 3 present additionally two more functions — cubic and linear, demonstrating similar described behavior. 
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S 
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Fig. 2. Comfort function dependence k(x) = 
1,x>1 
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Fig. 3. Comfort function dependence k(x) = 
1x>l1 


Indeed, suppose we need to embed an object of size / with comfortable distances on the left D_ and right D, into the 
space bounded by size L, and conditions /< L< D_+/+D, are met. Then, if one-sided comforts are presentable in the 
form K_(X_/D_)=k(x_), K,(X,/D,)=k(x,), where x_=X_/D_, x, =X,/D, , and dependence k(x) satisfies the 
above conditions, then, from the minimum condition of the objective function (2) we obtain: 

K_=K,=>k(X_/D_)=k(X,/D,). 

For the given nonlinear comfort function k(x), the resulting equation can be solved numerically by one of the known 

methods. However, since one-sided comforts are described by self-similar function k(x), from equation 


k(X_ / D_) = k(X, / D,) , we obtain relation X_/D =X,/D,, from which, taking into account equality 
L=xX_+1+X, , equations follow (3). 


Thus, it is shown that simple and rather convenient equalities (3), which provide embedding an object with optimal 
comfort, are valid for any one-sided comfort function k(x). 

Research Results. In a real situation, there is a need to place several objects at once. In this case, it makes no sense 
to solve sequentially the problems of optimal placement of the first object, then the second, third, etc., since when placing 
the next object, a need arises to shift previously installed objects so that the placement is optimal for the totality of all 


objects. If condition L< 310 is met, the objects in principle do not fit in the free space of length L. If 


i=l 


Le s ( DO +10 + D” ) then the objects can be placed so that they do not interfere with each other. In reality, the problem 
i=l 


of optimal placement arises if the following conditions are met: 


YIM <L< (oP +1 4p ). (4) 
i=l i=l 
In this case, we introduce the target function: 
1 2 n 
K,=KO+K4+..4K™, (5) 


where KY — comfort of the i-th object determined by formula (2). Here, when there are two consecutive objects with 
numbers 7 and i+1, then their neighborhood will be comfortable if the distance between them is not less than Dp! + Dl'*!) 


, which corresponds to the new comfortable distances dD = Dl) =D” + D('*) , i=1,2,...,(n-1) since a comfortable 


distance to the wall is one thing, and to another object, from where something can be pushed, is quite another. The rational 


arrangement of embedded objects is determined by the minimum of the target function (5). 
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The minimum of the target function (5), taking into account (2), gives us a system of n equations: 


K® =0; K® =0; ..., K& =0. (6) 
From (6), it follows: 


K =K; §=1,2,3,..50. . 


The system of equations (7) can be written in the following form: 


x x?) 
‘om “air 
(i) (i+1) 
ee eg |, eG, 
dD” B® 


(x) Ree ® 
| ae | — 2 I, 


That is, we have obtained n equations with respect ton unknowns X(), X),...,.X(") , where X(!) — distance of the first 
object from the left edge of the embedding area, X('), i = 2,3,...,2 — distance between objects with numbers i and (i-1). 


Since function k(x) is monotonic, system (8) is reduced to a linear system of equations, which does not depend on the 


type of the comfort function itself k(x). 


x9 x 
Uae aan, 
BO BO” 


Sees i=l 
(n) dD” . (9) 


- + 


System (9) can be solved by one of the known methods. However, due to the fact that the matrix of system (9) is 
highly sparse, the solution to the system can be found quite simply. For example, in the first (n-1) equations, it is possible 
to express X('#) by X( in each i-th equation, then, substituting this into the last equation, to obtain a linear equation 


with respect to X(). After that, moving from the first equation to (n-1)-th, we can consistently find X), X@),...,X™. 


In the case of placing two virtual objects, n=2, system (9) takes the following form: 


au x?) 
DO Be 
x9 pexO_x@ _79 _]72) 
52) A") (10) 
From system (10), we find: 
(z —1) = 1°) Bp” 
x) = eee 
"(B+ Dd?) + DBM 
11 
(L-19 1°) BS 1) ot 
x) = + 


(B° + 0°)" +d” bo 


where X w — distance between the first object and the left edge of space; X @) — distance between the right edge of the 


first object and the left edge of the second object. 
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Consider example (12) when: 


L=100, 1 =40, 1° =30, D® =10,D" =5,D° =11,D" =13. (12) 
Since conditions (13) are met, then: 
19) 419) <£<19 4194 DO +d) 4D? + po (13) 


For optimal arrangement of two objects, we can use expressions (11), which give regardless of the form k(x) in this 


case X") = 7.692; XO x 12.308; KY) = K?) =0. Here, the value of one-sided comfort depends on the type k(x). For 


linear dependence k(x), shown in Figure 3, K® =x" = K”) aK” ~ 0.769. If dependence k(x) corresponds to 


Figure |, we get the comfort value Kk” = KH = K”) = KO ~ 0.973 . Thus, the paper introduces the concept of optimal 


placement of a set of virtual objects in physical space. A model has been developed that provides solving the problem 
of optimal placement of a set of objects. It is shown that the solution to this problem does not depend on the type of 
monotone comfort function. 

Discussion and Conclusion. The theoretical aspect of the important issue of optimal placement of a set of virtual 
objects in physical space (the problem often encountered in augmented reality applications) was considered. By 
proposing a new mathematical model and including fuzzy logic, we laid the foundation for an algorithm that could 
potentially help users find a rational and convenient location of virtual objects in their real environment. 

The foundations laid in this study show that the proposed model effectively solves the problems associated with 
the placement of a set of virtual objects in a given physical space. By analyzing virtual planes and taking into 
account the distances between virtual objects and the edges of these virtual planes, our method provides optimal 
placement considering the linear dimensions of virtual objects and the comfort zone around them. 

The results obtained contribute to the current development of augmented and mixed reality applications, 
providing a theoretical solution to the problem of optimal placement, which, in turn, can improve user interaction 
and overall satisfaction with the tools of the technology under discussion. Moreover, the possible applications of 
this research go beyond AR applications, they open up new avenues for research in the related fields, such as virtual 
reality, mixed reality, and spatial computing. 

Considering the results of this study, future developments may be aimed at verifying the algorithm through 
empirical testing, enabling dynamic real-time adjustments based on user behavior, and exploring the integration of 
our approach into various XR application scenarios. As the area of augmented reality continues to evolve, we expect 
that our research will make a significant contribution to the development of the technology, inspiring its wide 


dissemination and further enriching the user experience. 
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Abstract 

Introduction. In the field of computational mathematics, there are many ways to approximate the model of fluid 
mechanics. Methods and estimates of approximation quality criteria, such as stability and convergence, are developed, 
while a combination of approaches to constructing economical difference schemes, such as splitting by physical 
processes, regularization by B. N. Chetverushkin, a linear combination of the Upwind and Standard Leapfrog 
difference schemes in aggregate has not been implemented and evaluated before. The authors were faced with the task 
of approximating each part of the hydrodynamic model split by physical processes with the most adequate scheme and 
further investigating the correctness of this approach. 

Materials and Methods. The mathematical model of hydrophysical processes is closed by the empirical equation of 
the state of salt water. Significant properties were selected, a mathematical model was built. Difference operators 
approximated differential operators. An algorithm for layer-by-layer modeling of transients was constructed. The 
algorithm has been implemented in the form of the program, which mainly contains elementwise (massively -parallel) 
operations. 

Results. Mathematical models of hydrodynamic processes in reservoirs were obtained, taking into account three 
equations of motion in the presence of a density gradient of the aqueous medium when hydrostatic approximation was 
abandoned. A new method of calculating the pressure field using B.N. Chetverushkin’s regularizers in the continuity 
equation was tested. A software module for numerical simulation of hydrophysical processes of water movement with 
different salinity and density was developed. This is open-source software that provides not only the redefinition of 
empirical dependences (as algebraic functions), but also the connection of external simulating modules to display 
dependences algorithmically. 

Discussion and Conclusion. The developed model of hydrophysics, taking into account the properties of salt water 
and the dynamic relationship of the mechanical movement of water with salinity, can be used to study the formation 
of a nonequilibrium distribution of parameters and identify the most stable parameters of the aquatic environment. The 
model explains the downward movement of oxygen. That will help in the future to estimate the values of the parameters 
of the aquatic environment, which are difficult to measure directly. It can be used in the procedure of parametric 


identification of hard-to-measure parameters of the aquatic environment. 


Keywords: mathematical model, stratification, seawater dynamics model, quasi-hydrodynamic model, cell occupancy 


method, central difference scheme, sweep method, FTCS scheme, Upwind Leapfrog, Standard Leapfrog 
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AHHOTAalnaA 

Beedenue. B oOnacTu BEIYHCIUTeIbHOM MaTeMaTHKH VM3BeECTHO MHO2KECTBO CHOCOOOB alllpoKcHMalMu Moje 
MeXaHHKH %KH[KOCTH. YUeHbIMH BbIPAOOTAHbI MeTOAbI HM OWCHKH KpHTepveB KayecTBa AalINIpOKCHMallHu, TaKHX Kak 
yCTOWYMBOCTb HU cxogqMMocTb. KoMOuHalluaA MOAXOLOB MOCTpOeHHA Z9KOHOMHMYHBIX Pa3HOCTHBIX CXe€M, TaKMX KaK 
pacijemeHnve mo du3vyeckKHM mpoleccam, perysapu3sauua mo b.H. Yersepymikuuy, sMHetHad KOMOMHaMA 
pa3HOCTHOM CXeMbI «KaOape> UH «KpeCT» B COBOKYIHOCTH paHee He peasIM30BbIBalacb MU He OWeHUBalacb. Ilepey 
aBTOpaMH CTOAa 3aaua alllpOKCHMMpOBaTb KaxKAy!O YaACTb paculeréHHOM MO (u3H4eCKHM Mpoweccam Moe 
THpOAHHaMUKU HanOosee aleKBaTHOH CXeMOM HU aylee MCCICOBaTb KOPpeKTHOCTh JaHHOro MOAXoMa. 
Mamepuanvi u memoodvi. Matematuueckast MOJeIb THApOPu3sH4eCKHX MpOLeccoB 3aMbIKaeTCA IMIMpH4ecKHM 
ypaBHeHHeM COCTOAHHA comeHow BOAbI. BrlOuparoTca 3HadMMbIe CBOMCTBa, CTPOHMTCA MaTeMaTHYeCKad MOJEIIb. 
Pa3HOCTHbIe OMepaTOpbl alllpoKcHMUMpyloT AupmepenuMalbHble oMepaTopsl. CTpovtca alropuTM mocsoMHoro 
MOJ{eIMpoBaHHA TepexoJHbIxX Upoueccos. ANTOpHTM peasIM30BaH B BHe MpOrpaMMBbl, KOTOpasd, B OCHOBHOM, 
COJepXUT MOSIEMeHTHbIe (MaCCHBHO MapasWIesIbHbIe) OllepalHu. 

Pesynomamoi uccaedoeanua. Wlonyaenbl MaTeMaTH4eckHe MOJeIH THAPpOAMHAMMYeCKHX TIPOWeccoB B BOJOeMAaX, 
YUHTEIBAIOMHe TPH ypaBHeHHA ABMKCHHA IPH HaMYMU TpayMeHTa MIOTHOCTH BOAHOM cpeAbl Mp OTKa3e OT 
rHypoctaTuyeckoro upHOmMKeHHa. AtmpoOvpoBaH HOBbIM ciocoO BBIYMCIeHHA MONA WaBeHua C MpHMeHeHHemM 
perynapu3atopos no 5b. H. YetsepyuikuHy B ypaBHeHHuH Hepa3ppiBHocTH. Pa3padoTaH porpaMMHbIii MOJlYJIb 
YMCICHHOTO MOJeMpOBaHHA THApOPU3M4eCKUX IIpOWeccoB JBHXKCHHA BOAbI C pa3IMYHOM CONEHOCTBIO HU 
IVIOTHOCTBIO. ITO OTKPbITOe IporpaMMHoe ObecrieyeHHe, AOMycKarollee He TOKO MepeonpeyeueHve IMIMpPHYeCKHX 
3aBHCHMOCTeH (Kak allreOpamueckux dyHKUMi), HO HM MOAKIOYeHHe BHeEWHUX MOJeIMpyIOWMX MOAyIeH DA 
oTOOpaxkKeHHA 3aBMCHMOCTeH aIrOpHTMHYeCKH. 

O6cyacdenue u 3akjiio4uenue. Pa3padoTaHHad MOJeb THApOPUu3HKH, YYMTHIBAaIOWad CBOMCTBa CONEHOM BOJbI H 
JMHaMHYeCKYIO CBA3b MCXaHHYeCKOTO ABYWKCHHA BOLI C CONGHOCTbIO, MOXKET MPHMCHATLCA WIA U3yYeCHHA 
(bopMvpoBaHHA HepaBHOBeCHOro paciipeseeHHA MapaMeTpoB HU UACHTUPUKauHU HanOoee CTAOMIbHBIX NapaMeTpoB 
BOAHONM cpegbl. Moyen oObACHAeT HUCXOAIee WBYXKeHHe KHCIOpOza, YTO NO3BONHT B OYyAYUIeM OLeHHBaTb 
BeJIMYHHbI WapaMeTpOB BOHOM cCpebl, KOTOPbIe CJIOX%KHO U3MeEPHTb HelocpeyCTBeHHO. OHa MOxKeT OBIT 
MCIOJIb30BaHa B Ipolleype NapaMeTpH4eckon HACHTHPUKALMH Tpy {HOU3MepsAeMBIX MapaMeTPOB BOAHOH cpesBl. 
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Introduction. One of the critical tasks related to the ecology and life safety of people living in coastal areas is 
forecasting and modeling the movement of water in the seas and large regional reservoirs. Moreover, it is important to 
study the process of transport of substances dissolved in an aqueous medium, taking into account stratification and the 
dependence of water density on many variable factors. Such predictive modeling is likely to make it possible not only 
to assess water quality, but also to manage it under conditions of climate change and industrial impacts. A more general 
and major goal is to link water quality to the number and species diversity of natural hydrobiocenoses living in the 
hydrosphere. When solving these problems, it is necessary to take into account the hydrodynamic characteristics of 
the aquatic environment, features of external influencing factors, such as the heterogeneity of temperature distribution, 
salinity, oxygen saturation of water, the amount of gases dissolved in water, acidity. These parameters are some key 
parameters of the biological activity of the aquatic ecosystem [1]. 

The water area of the reservoir can also be considered as a transport system that transfers oxygen from the 
atmosphere to bottom sediments. However, there are cases of the formation of chemical gradients of large magnitude 
at relatively small depth differences — in boundary layers with a relatively small value outside these zones. The reason 
for the appearance of such sites with the general stratification of water by density in universal gravity on the one hand, 
and the radiation effect of the sun on the water, causing its heating, on the other. These processes can cause a decrease 
in the rate of production, destruction and recycling of biogenic elements and bio-organisms until these processes stop, 
as well as predetermine the biodiversity of hydrobionts in general and species composition in particular [2]. 
Temperature stratification affects significantly the distribution of organisms in the water column, the transfer and 
deposition of impurities harmful to bio-organisms. An increase in the temperature of surface waters causes a violation 
of vertical water exchange and, accordingly, a decrease in aeration of the deep-water zone, a decrease in solubility and 
oxygen concentration in water. Stratification by density, temperature and chemical composition limits the convective 
rise of biogenic elements, carbon dioxide and products of incomplete oxidation of organic substances entering the 
hypolimnion (cold, salty, dense layers of water) as a result of sedimentation of seston into the surface layers of water. 
From the beginning of stratification until its end, the surface layer is depleted, and the hypolimnion, on the contrary, 
is enriched with these substances. As a result, physico-chemical stratification causes an uneven distribution of a number 
of biologically significant substances in depth and is the reason for the self-organization of a complex structure of 
ecological niches [3]. Mathematical modeling of mechanical, chemical and biological processes occurring in aquatic 
ecosystems is urgent, it is associated with problems of ecology and life safety of the population of coastal territories. 

We should name the outstanding scientists who made a significant contribution to the study of hydrology and 
oceanology. V.P. Dymnikov was engaged in the study of climate and oceanology, modeling of the atmosphere and 
ocean. A.S. Monin and M.Y. Belevich investigated the processes of kinematics of the aquatic environment, turbulence 
and microstructure of the ocean. Soviet scientist A.S. Sudolsky studied the dynamics of waters and coastal processes 
in various reservoirs in relation to solving the problems of designing, building and operating specific hydraulic 
structures, and rational economic use of reservoirs in general [4]. According to V.I. Vernadsky, one of the most 
important manifestations of life is the gas exchange of organisms and the environment, mainly respiration processes 
based on oxygen consumption [5]. Some outstanding Russian scientists, including S.V. Bruevich, whose works were 
devoted to the development of analytical research methods, the formulation of the basics of hydro- and bio- 
hydrochemistry, were engaged in the study of hydrobiological processes of reservoirs. G.G. Matishov and V.G. 
Ilyichev actively study the conditions of optimal exploitation of water resources, develop models of transport of 
pollutants in water bodies, and investigate the assessment of their impact on the bioresources of the aquatic 
environment [6]. 

To study the influence of these processes of the aquatic environment, a complex of interrelated mathematical 
models is being developed. It is based on the use of accurate predictive models and software implementation of 
economical numerical methods, which provide detailed investigation of the kinematics of the process, cause -effect 
relationships, and the state of the modeling object. The existing methods and means of predictive modeling of the state 


of the aquatic environment, taking into account a number of biotic and abiotic factors, including the processes of 
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distribution of oxygen, carbon dioxide, salts, are based on general scientific approaches, simplified mathematical 
models with low adaptability, lack of modeling of nonlinear dynamic processes characteristic of most aquatic 
ecosystems, incorrect setting of the boundaries of the computational domain. In some cases, the studies are accompanied by 
a formal definition of boundary conditions, and give rather rough and approximate simulation results. 

When modeling the substance transfer process based on the advection-diffusion equations, a good approximation 
of advective terms, which are gradients of pressure, mass density and total energy, momentum of motion, is required. 
The use of standard difference schemes with large estimates of similarity parameters causes loss of calculation accuracy 
due to an increase in the approximation error and increased restrictions on the time step due to the stability condition 
of the difference scheme. In the works of A.I. Suhinov, A.E. Chistyakov and others [7, 8], itis shown how to effectively 
use a linear combination of the Upwind and Standard Leapfrog difference schemes with the optimal values of weight 
coefficients for the approximation of the transport equation. The efficiency of these methods is achieved by 
optimization of the approximation error of discrete-continuum model, the exact solution of the mass transfer with 
constant speed. Studies have shown that this approach also extends to hydrodynamic models with a variable (sign- 
alternating) velocity without the effect of grid viscosity. Another positive quality of such approximations is that they 
can be used to simulate complex flow structure, e.g., vortex. Currently, numerous researchers use such a scheme for 
the simulation of turbulent flows. Members of leading foreign research organizations, such as Stanford University, 
Imperial College London, etc., as well as members of the Institute of Computational Mathematics of the Russian 
Academy of Sciences E.M. Volodin, A.V. Glazunov, A.S. Gritsun, N.G. Yakovlev and others [9] published works in 
which mathematical modeling of climatic changes, hydrodynamic and atmospheric processes and phenomena are 
carried out on the basis of eddy-resolving schemes. The quasi-hydrodynamic approximation of a continuous medium 
enables to further reduce the time step requirement and increase the spatial resolution of the model with limited 
computer memory. In practice, a small term proportional to the second derivative in time or density is added to the 
system of Navier-Stokes and continuity equations. This approach makes it possible to smooth out non-physical 
fluctuations in mass density and momentum, as well as total energy, brought along the spatial grid faster than the speed 
of sound. 

Existing universal application software packages (e.g., “Mars3d” software package, Ecointegrator, CHARISMA 
software package, SALMO complex, CHTDM, CARDINAL software package, packages for modeling various 
processes of aerohydrodynamics, PHOENICS, FLUENT, GAS DYNAMICS TOOL software packages) do not take 
into account some properties of the simulated complex systems, thus reducing the accuracy and efficiency of modeling. 
These properties include spatial heterogeneity of the aquatic environment motion, vortex structures of currents. 
Mathematical models and algorithms for their numerical implementation do not take into account the probability of a 
significant change in the depth and density of the aquatic environment, which can cause instability of the numerical 
solutions obtained. For this reason, such specialized software packages can be used to simulate a limited variety of 
hydrophysical processes of water systems. Most of the well-known specialized software (ADAM, CAL3QHC, Chensi, 
TASCflow, ISC—3, PANACHE, REMSAD, UAM-IV, ECOLOGIST, PRISM, VITECON), designed to simulate the 
process of spreading pollutants, the interaction of hydrobionts, is mainly focused on uniprocessors, represented mainly 
by personal computers. In such systems, only single component modules of these systems (e.g., ECOSIM and 
MAQSIP) are scaled to parallel systems. In practice, a small term proportional to the second derivative in time or 
density is added to the system of Navier-Stokes equations and continuity. This approach makes it possible to smooth 
out non-physical fluctuations in mass density, momentum and total energy brought along the spatial grid faster than 
the speed of sound. 

Materials and Methods. The success of the development of a mathematical model of hydrophysical and 
hydrobiological processes depends on the availability and elaboration of test cases and tasks for investigating steadily 
observed phenomena in the seas, such as vertical mixing and redistribution of salinity and oxygen, halocline and 
thermocline. To study these phenomena, the paper uses a model of hydrodynamics that takes into account the balance 
of mass forces and cross-border flows [1, 2]: 

0 
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k+ 
oe ce div(Tv)+(v,pb)+(V.h)+pq, 
ot. p 
where pv — flow density; Ok / Ot — kinetic-energy density rate; 0s / Ot —-— self-energy density rate; T — stress 


tensor Tn =—pn; b — mass force; h= A(t, *) — heat flux density; g — specific heat inflow due to radiation. Due 


to the fact that these phenomena are usually described by a vertical distribution of parameters over depth, it is 
reasonable to obtain a simplified model that provides rapid identification of unobservable parameters. We select a 
cylindrical region V with bases on the bottom and surface of the water with a cross-sectional area S. We project the 
velocities, flows and forces to the vertical direction, assuming that the partial derivatives in x, y of the parameters to 
the horizontal direction are equal to zero outside the cylinder, i.e., assume horizontal uniformity of the parameters of 


the aqueous medium. Let us write system (1) in a conservative form so that it enables to determine the mass density 


(p), mechanical momentum (pv ) and total energy density (p(k +8) ), k=v?/2. In the cylinder, we allocate 


infinitesimal volume and assume that each such volume, which is V, is affected by the reaction force of the support 
(bottom of the reservoir), equal to hydrostatic pressure, and a force similar to the friction force caused by the viscosity 
of the fluid and momentum transfer, not equal to zero with vertical movements of the fluid. 

If we neglect horizontal fluid movements and assume that only vertical movements are essential, and the fact that 
the density of the aqueous medium significantly depends on salinity, accept the simplifications and conventions 
outlined earlier, then the equations of hydromechanics [1—3] in compact form are presented by a system of partial 
differential equations [2]: 

ap , (or) 
Ot Ox 


f) 


A(pv)  Olpv?) __ a 
Ot Ox Ox 


of a 
+F/S,F =—pgS C )e=f0.0 
Ox\ Ox 


a ; a (pv) Fv. 6(, aT 
aL (ety /2)]+=[v(p(e+v?/2))] = ae rok =| k(T -T,,,),T = f(e)=e/c, (1) 
as (vs) 2 “| 
ae = uu 5 
ot Ox Ox ox 


where p= f(p,7,5), €=f(p,v,T) — the empirical equation of the state of seawater and the equation, closing the 
system by internal energy, respectively; p — total hydrodynamic pressure; p — local density of the aquatic 


environment; v — projection of the velocity vector function on the vertical (z axis is directed upwards from the bottom 


to the surface); ¢ — volume density of internal energy; p — pressure of the gas enclosed in an elementary volume 


between adjacent layers; F — volume density of the generalized force (sum of forces) applied to elementary volumes 


of fluid, in addition to pressure; s = s(x,t) — salt concentration; S — cross-sectional area of the cylindrical selected 


area in which the most intensive upwelling and salt transport process is assumed; 7 — absolute water temperature; 


T,,, —-temperature of the water external to the allocated volume; k — thermal conductivity of water; 


env 
g =9,8m/ s? — acceleration of gravity; 1. — coefficient characterizing the intensity of momentum transfer due to 
viscosity. 
The boundary conditions characterizing the property of impermeability of fluid through the rock that makes up the 


bottom of a water body for system (1) can be written in the form of equalities: 


ao a 5, SF 0. (2) 
on on on 


Where n — the normal vector directed inside the computational domain. 
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The International Thermodynamic Equation of Seawater [10] defines density p of seawater as a function of salinity 


Ss, temperature T and hydrostatic pressure p, which has the form: 


p(s,T,0) 


1-p/K(s,T, p) @) 


p(s,T, p) = 


where K(s, T, p) — average modulus of elasticity; digit 0 corresponds to one standard atmosphere (101,325 Pa). 


Under numerical modeling based on (1) — (3), also in demand as (3), we establish an algebraic relationship of 


pressure with density, temperature and salinity. The quadratic equation has two roots: 


_ VD +(1-A)-p(s.T, p)+ A: py (8.7, p) 


Pir D, : (4) 
Dy = 2B-(p(s,T, p) —Po (sf, P)) 
D=(-4BK,(s,T, p)+ A? —2A+1)-p(s,T, p)? +(8B-K,(s,T, p)—2A? +2A)- p,(s,T, p)- (5) 


p(s,T, p)+(A* —4B- Ky (s,T, p))- py(s,T, p)?; 


where A, B, p, (s,T, P), K, (ss T, P) — variable parameters whose relation to temperature, salinity and density is 
determined by the standard model of seawater; p(s, T, p) =p, where T — water Celsius temperature. The speed of 


sound in salt water can be expressed by the formula: c(s, T, p) =Op/op. 


Suppose that the momentum, salinity and heat of a volume of water, bounded by a cylindrical surface, change only 
when the parameters at the boundary of the domain (Dirichlet boundary conditions) change, in particular, the salinity 
changes on the surface of the water. The behavior of mechanical parameters — density, momentum — is set by 
Neumann conditions and is inaccessible to external influence. Regardless of these assumptions, within the boundary 
of the aqueous medium, the system of equations (1) can be written in the vector form of a transfer-reaction model with 


a source additive [11, 5]: 


—+—=T, (6) 


env 


where T = (0 F/S+ty% Fv/S+k0°T/dx? —k(T -T. yy — vector characterizing the interaction of the flow with the 
surrounding fluid and the planet; gq=(qgpvE) — _ vector of conservative state variables; 
f= (pv pt+pv? v(E + p)) — vector of flows acting as feedback and closing the balance equation. From the 


components of vector f common multiplier can be derived f =v-q. 


In practice, to solve the transfer equation of form (6), a difference scheme has proven itself well, whose layered 
transition operator is obtained by a linear combination of similar transition operators of the Upwind and Standard 
Leapfrog schemes [12—15]. Taking into account the relationship between mass flows and density changes turned out 
to be more effective for a quasi-hydrodynamic system regularized according to B.N. Chetverushkin, approximated 
according to the FTCS scheme (forward-time central-space). 


The discretization of the continuum model will be carried out by the integro-interpolation method on a uniform 


grid S=S,xS, S, =1%, =ih,i=0+n,n-h=L} , 55 = {ti = jt, j=0+m,m-t=T} , where h — vertical grid step, 


i — index of the node (control/final volume in terms of the Godunov method) when numbering in space, t — grid 
step in time, j — the number of the time layer. We describe a finite-difference approximation of the model. The 


difference schemes used to approximate the first-order balance equations give relatively acceptable results only with 
a very small grid step [3], which causes an active consumption of resources of computing devices and developers who 
create algorithms and programs that numerically implement them. 

The problem of setting boundary conditions and their coordinated assignment at alternating speeds was solved by 
the method of filling control cells [7]. Taking into account the maximum occupancy and expressing its functional 


dependence on the node number makes it possible to increase the accuracy of approximation of boundary conditions. 
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Investigations on the theoretical and experimental selection of approximation methods and schemes led to the 
choice of splitting methods [13] and regularization. The continuity equations were regularized by B. N. Chetverushkin. 
The given equation was subjected to a difference approximation according to the FTCS scheme. The balance (transfer) 
equations of momentum, salinity and total energy were replaced by explicit equations obtained by linear combination 
of various approximations of the transfer operator (leapfrog and Upwind Leapfrog) [12, 15]. According to the splitting 
method, the transfer of momentum (and velocity), salinity, and total energy was approximated in a fractional step. The 
water velocity variation, which determined the intensity of transport and mass flows, was introduced at the second 
fractional step by solving the wave equation. Denote p=p,, p=p"*°, p=p™!,vev,, =v, Pav", O<oK<l. 
In the split form [16], after the introduction of a regularizer into the continuity equation and semi-discretization on 


grid S,, approximation of the partial derivative in time by the FTCS scheme, the first two equations of system (1) can 


be written as: 


py—pv O{ ov 
PYERY aise agen : 
o SUA DE of =), (7) 
1 p= p: mer! Sil Pett —2 Pa + Prt _ 7 
-( 7 +(68), Px Cr } f(p.T),c=c(s,T.p), (8) 
ne (9) 
T 


Equations (7) and (9) are a discrete analogue of the momentum transfer and change model — the second equation 
of system (1). Equation (8) predicts a change in the pressure field taking into account the continuity of the flow and 
information about internodal (intercellular) and boundary flows. Assuming time constancy at each step of the speed 
of sound, equation (8) can be approximated by an implicit time difference scheme and solved by sweeping. It was this 
direction in numerical modeling that was selected for the study. 

Equation (7) and the fourth equation (1) are approximated on an equidistant grid by an explicit difference 
scheme [13]. Summands (8) determining the wave properties during the propagation of the pulse (the right side of the 
equation) are approximated by the central differences in space: 

kaj prt —2pt + pr! ( , Phicee | pitt pe 
C2 2 Ii h2 ai h2 ; 


(10) 


where k, — occupancy of the cell to the right of the node with index i; k, — occupancy of the cell to the left of the 
node with index i; k, — degree of occupancy of the control area X,_,,. SX S$ X;,;,., located in the neighborhood of the 
node with index i(k, ,k, ,k, — iindex-dependent variables) [7]. 

The remaining two summands (8) determine the relationship of the rate of change in density with its flow. (8) includes 
the ratio of the rate of change in density to the time interval Tt, expressed by the difference with a forward shift: 


Ti 
PP -P? 
T 


k 


0.1 


the mass flow is transferred to the right side (8) and approximated by the central differences: 


1 0(pi sea staill: =. 
ky, 1 (pi) ~k,, ((PU); 1. —(PU);) +k,, ((pu), — (pu); 1.) (ib 
"Tt Ox | h-t h-t 
where k,,k,,k, — degree of occupancy of the areas located in the vicinity of the cells with number i[7]; 
k, characterizes the occupancy of the area [x,_,,%,,,]; k, — [%,,%,], k&, — [%.,.%]. 
Balance equations of the mechanical impulse and pressure (9) are approximated by FTCS: 
yet! _ yet n+l _ n+l n+l _ n+l 
k,,? i p i = k,, Pe! Pj ee Pi ’ (12) 
, Tt 2h 2h 


Let us present the problem of solving the difference approximation of equation (8) in the form of a matrix sweep 


problem with a vector variable in time on the right side Ax=F: 
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AX, +C,x, + BX. = F,. 
The approximation of equation (8), defined on a three-point difference template, in the form of a linear system of 


equations Ap"! = f(p", p",(pv).) , Solvable with respect to pressure (p), has the form: 
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A similar band matrix is obtained during a similar approximation of the other two terms (8). 

The model can take into account dynamic changes in the flow at the boundaries and other parameters, such as 
temperature, salinity, oxygen content. Therefore, the simulation of water movement must be performed in a time cycle 
in layers with the retention of operational information about the parameters on at least two time-adjacent layers of the 
solution to the grid equation. The algorithm for calculating hydrodynamic parameters on a two-index grid in space and 
time includes: 

— building of the projected pulse changing according to the first equation (3); 

—approximate calculation of the function of the spatial distribution of pressure as a function of density and 
temperature from the previous time layer; 

— assessment of the density changes according to the second equation (3); 

—calculation of the pressure gradient from the new values of density and temperature, correction of the momentum 
distribution according to the third equation (3); 

—finding a new distribution of the total energy and temperature according to the approximated third equation of 
the system (1). 

When algorithmizing the method for solving a system of split equations, the following are introduced: 

—binary-numeric masks that predetermine the switching of the difference scheme template with a change in the 
sign of the velocity; 

—variable shifts of the indices of neighboring nodes, which make it possible to write the solution to equation (12) 
for boundary and internal nodes by a single system of computational operations. 

Such variable index shifts are used in calculating gradient approximations both in the volume of the continuum 
model and at the boundaries with the second order of accuracy according to discrete analogues of equations (7) — (9) 


approximated by FTCS and by a linear combination of Leapfrog and Upwind Leapfrog schemes [12-14]: 


n+l _ qn 4 no _ qn-l n—qn n—gqn-l n 4 qn 
qi qi se qi-1 Gi fs ve qi qi-1 a qi qi ef vn Givi Gi-1 _ 0, vn > 0, 
Tt 3 2t h 3t 3h 
ntl — qn 4 n _ qn-l n an n —qn-l nao _ogn (14) 
qi Gi) Gist Gin yn Gin 7G |, GWG yn Gin ~ Gia _ O,v" <0, 
v 3 24 h 3t 3h 


where is characteristic q for transfer of salt (q =s ) and momentum g=pu. 


Parameter m is a “switch” of the flow when the sign of the velocity changes in the difference approximations of 


the pressure gradients and the mass flow: 
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0, v<0, 


fi=1—-mm={? v>0. 


When using such a switch, the approximation of the gradient in solving the transfer problem for momentum and 


salinity can be written by the formula: 


0q 
Ox|,_ 


w2 mf_#a am fn 4 4h Gis ~ Fins ; (15) 
3 h h 3. 2h 


The use of the cell occupancy method makes it possible to correctly approximate the boundary conditions when 
placing water state parameters in computer memory with an array of numerical values, and variable index shifts enable 
to reduce the volume of symbolic recording of the subroutine for performing layer-by-layer iterative changes of state 
variables assigned to the nodes of the difference scheme [17]. 

The matrix equation with a band tridiagonal matrix A (13) is solved by sweeping [18]. For the initial verification 
of the software implementation, a simplified model of the state of salt water in the form of P. S. Lineikin’s equation 
was used. 

Research Results. The simulation was performed on the basis of a software package written in MATLAB. 
Debugging was performed using the GNU Octave interpreter. Operation also presupposes the availability of libraries 
of this system. The software package consists of 16 functional modules. The selection of this interpreter and the 


corresponding language was due to the ability to write and test a program that operates with an array of state variables, 


which is a projection of the desired functions p(x,t) , v(x,t), &(x,f) on a spatial grid. At this stage of the study, the 


authors abstracted from the specifics of performing element-by-element operations on arrays of real numbers. 
The software system consists of the interconnected subprograms: 

— execution of one step on the transfer equation according to formulas (10); 

— calculation of seawater density (unesco_urs); 

— calculation of the pressure of seawater located in the gravity field at a given temperature, salinity and density 
according to empirical equations of state (rhoTS2P); 

— calculation of the speed of sound depending on temperature, salinity and pressure according to the standard and 
the updated UNESCO formula (speed_of_sound); 

— estimation of the fluid viscosity when its volumes move under the influence of pressure forces at all points of the 
spatial grid (ForceOfFriction); 

— cyclic variation of state and time variables during vertical movement of salt water (aqua_process); 

— initial setting of constant values characterizing the fluid, initial and boundary conditions of its global 
hydrophysical equilibrium (start, set_parameters); 

— formation of band matrix approximating the equation containing pressure (func5); 

— solution to the matrix equation by the sweep method (run_sweep_shuttle); 

— temperature <> total energy conversion (TFromE, EFromT); 

— calculation of total energy balance (TotalPower); 

— solution to the diffusion-convection-reaction problem, including the approximation of a combination of Leapfrog 
and Upwind Leapfrog difference schemes (ADR_solver); 

— estimation of gradients and time derivatives with respect to central and directional differences, taking into account 
the change in the sign of the velocity and the pattern of the difference scheme (diff123). 

The constructed model was used in a test run to estimate the density change under an increase in the salinity of the 
surface of the aquatic environment by 1% at the 30th second from the start of the simulation and a maximum depth of 


1 km. Density transients in a set of cross-sections (function p(t, x)| ;)» Spaced from each other, caused by an 
xeS, 


instantaneous change in salt concentration, are finite in time (Fig. 1), and the density pulse front is greatly weakened 


when moving in space. 
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Fig. 1. Graphs of dependence of density on time with a sharp change in salinity on the water surface 


Discussion and Conclusion. The presented spatially distributed model does not yet allow predicting a steady 
redistribution of density and a change in the salinity gradient in numerical calculations, because it does not take into 
account the buoyancy of less salty water and its change during salinization of the upper layers. Mathematics and 
software, which greatly simplify the modeling of processes leading to the observed effects of halocline, chemocline 
and pycnocline, has been tested and debugged on a variety of test tasks (momentum and salt transport, pressure wave 
propagation). For high-precision modeling of the observed physico-chemical phenomena in the seas, it is required to 
solve additional tasks of identifying model parameters, taking the data of observations and remote sensing of the Earth 


as initial information. 
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